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Streszczenie w jezyku polskim 


Tematem niniejszej pracy jest analiza probabilistyczna wrażliwości mechanicznej 
odpowiedzi konstrukcji inżynierskich na imperfekcje geometryczne oraz zmianę 
własności materiałów. Z tematem tym ściśle związany jest problem estymacji 
niezawodności konstrukcji, również szeroko opisany w niniejszej rozprawie. 

W pracy dokonano przeglądu metod wykorzystywanych w analizie 
probabilistycznej oraz zaproponowano procedury wykorzystujące te metody w analizie 
wrażliwości niezawodnościowej. Opisane zostały kolejno: metoda Monte Carlo (MC), 
metody redukcyjne próbkowania warstwowego (SS) i hipersześcianu łacińskiego (LHS), 
metoda estymacji punktowej (PEM) oraz metoda powierzchni odpowiedzi (RSM). 
Najwięcej uwagi poświęcono technice próbkowania ukierunkowanego (TRS — Targeted 
Random Sampling), na podstawie której stworzono program komputerowy w środowisku 
MATLAB, umożliwiający poszukiwanie próbek w otoczeniu stanu granicznego 
konstrukcji. Osobny rozdział poświęcono klasyfikacji i omówieniu technik analizy 
wrażliwości, m.in. analizy jednoczynnikowej, projektowania przesiewowego oraz 
analizy wariancji. 

Przedstawione metody probabilistyczne wykorzystano w analizie kilku 
przykładów modeli prętowych: kratownicy von Misesa, kolumnie Zieglera, wieży 
telekomunikacyjnej i kopuły prętowej. Wskazano optymalne sposoby określania 
wrażliwości mechanicznej odpowiedzi układu na zmiany poszczególnych zmiennych 
losowych. 

We wnioskach podsumowano przeprowadzone obliczenia oraz określono 
kierunki dalszych prac, których nadrzędnym celem powinna być implementacja metod 
probabilistycznych w standardowym projektowaniu konstrukcji. 


Streszczenie w języku angielskim 


This dissertation concerns a probabilistic sensitivity analysis of mechanical 
response of engineering structures to geometrical imperfections as and to material 
parameter variation. The issue is closely related to the problem of reliability-based 
assessment of structures, broadly addressed in this dissertation too. 

The thesis reviews the methods of probabilistic analysis and proposes procedures 
to apply them in the field of reliability sensitivity analysis. The thesis covers the following 
methods: Monte Carlo (MC), Stratified Sampling (SS), Latin Hypercube Sampling 
(LHS), the Point Estimate Method (PEM) and the Response Surface Method (RSM). 
However, most attention is focused on the Targeted Random Sampling technique, which 
constituted the basis for a computer program in the Matlab software environment, 
dedicated to samples search in the limit state vicinity of a given structure. A separate 
chapter deals with the classification and discussion of sensitivity analysis techniques, e.g. 
one at a time analysis, screening design and the analysis of variance. 

The probabilistic methods presented in this paper were illustrated by selected 
worked examples of bar structures: von Mises truss, Ziegler’s column, 
telecommunications tower and reticulated shell. Due to the analysed examples, alternative 
methods were pointed out to determine the sensitivity of the response to particular random 
variations. 

The concluding chapter summarizes different types of analysis and indicates 
possible directions for further work, directed to implementation of probabilistic methods 
in standard engineering calculations should be an overarching objective. 
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Rozdział 1. Przedmiot, cel i zakres pracy 


1.1 Przedmiot rozważań 


Proces projektowania oraz budowania konstrukcji inżynierskich powinien być 
przeprowadzony w sposób gwarantujący jej bezpieczeństwo oraz odpowiedni komfort 
użytkowania. Trzema podstawowymi warunkami, które musi spełnić każdy obiekt są 
kryteria: nośności, sztywności oraz stateczności. Innymi słowy, posługując się definicją 
zawartą w normie PN-ISO-2394, konstrukcja musi być zaprojektowana w sposób 
zapewniający jej niezawodność czyli zdolność do wypełniania funkcji przewidzianych w 
projekcie w danych warunkach eksploatacji Obliczenia deterministyczne, dominujące w 
inżynierskiej analizie konstrukcji, nie dają zadowalającej odpowiedzi dotyczącej 
niezawodności danego układu. Jedynie metody probabilistyczne, będące przedmiotem 
niniejszej pracy, pozwalają na uwzględnienie faktu podlegania przez wielkości 
wykorzystywane w projektowaniu wpływom czynników losowych oraz umożliwiają 
oszacowanie ich bezpieczeństwa. 

Dostępne i powszechnie stosowane inżynierskie narzędzia obliczeniowe nie 
pozwalają na przeprowadzenie pełnej analizy niezawodnościowej. Należy podkreślić, że 
w przeciwieństwie do programów inżynierskich wykorzystujących liniową a nawet 
nieliniową analizę konstrukcji, programy umożliwiające przeprowadzenie obliczeń 
probabilistycznych wymagają od użytkownika rozległej wiedzy z zakresu rachunku 
prawdopodobieństwa oraz statystyki matematycznej, przez co prawdopodobnie jeszcze 
długo nie zostaną w pełni zautomatyzowane. Na chwilę obecną inżynier nie dysponuje 
„czarną skrzynką” za pomocą której, po wprowadzeniu według odpowiedniego wzorca 
danych, uzyska wyniki bezpośrednio wykorzystywane w projektowaniu. Z uwagi na 
błyskawiczny rozwój możliwości obliczeniowych komputerów możliwe (a nawet 
konieczne) jest jednak włączenie metod probabilistycznych do standardowych analiz 
inżynierskich. Pytaniem otwartym pozostaje w jakim zakresie takie obliczenia powinny 
być wykonywane? 

Oczywistą odpowiedzią na to pytanie jest: określenie prawdopodobieństwa 
przekroczenia stanów granicznych konstrukcji lub utraty przez nią stateczności. Badając 
dostępne metody i algorytmy z łatwością można jednak stwierdzić, że właśnie takie 


obliczenia są najbardziej skomplikowane gdyż dotyczą tzw. ogonów rozkładów 
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prawdopodobieństwa, matematycznie trudnych do zdefiniowania i wyznaczenia. Nawet 
pobieżny przegląd literatury z zakresu niezawodności uzmysławia jak wiele różnych 
metod zostało zaproponowanych przy braku wiodących, powszechnie akceptowanych i 
uniwersalnych algorytmów. Właśnie ta mnogość możliwych rozwiązań świadczy o 
trudnościach obliczeniowych. Ponadto z uwagi na brak precyzyjnych danych 
wejściowych dotyczących rozrzutów wymiarów konstrukcji, parametrów materiałowych, 
a także (a może przede wszystkim) obciążeń, obliczenia niezawodności lub 
prawdopodobieństwa zniszczenia konstrukcji obarczone są dużą niepewnością. 
Zastosowanie metod probabilistycznych obejmuje szerokie spektrum zagadnień. 
Wśród nich istotne miejsce zajmuje analiza badań doświadczalnych materiałów i 
elementów konstrukcyjnych. Podejście takie umożliwia ich wszechstronny opis, obniża 
koszty badań, a także pozwala na uogólnienie opisu równań konstytutywnych lub 
mechanicznego zachowania się elementów konstrukcyjnych. Innym obszarem 
wykorzystania metod probabilistycznych jest zastosowanie budowanych obecnie 
algorytmów sprawdzających niezawodność istniejących, wzniesionych wiele lat temu 
konstrukcji. Dotyczy to szczególnie będących w nieustannym użytkowaniu obiektów 
zabytkowych, w których do czynienia mamy zarówno z nieuchronnym procesem 
degradacji wykorzystanych materiałów jak również ewolucją działających obciążeń. 
Należy podkreślić, że realistyczne ujęcie w obliczeniach procesów starzenia i korozji jest 
niezwykle trudne i podejście probabilistyczne jest w takich przypadkach szczególnie 
uzasadnione. Także rozwój produkcji nowych materiałów, np. kompozytów, z uwagi na 
ich niejednorodne cechy, wskazuje na konieczność stosowania opisu probabilistycznego. 
Kolejnym istotnym obszarem zastosowań metod probabilistycznych jest 
optymalizacja konstrukcji z uwzględnieniem jej niezawodności. Włączenie opisu 
niezawodnościowego do algorytmów optymalizacyjnych pozwala na ulepszenie procesu 
ekonomicznego projektowania konstrukcji z jednoczesnym uwzględnieniem 
niezbędnego poziomu jej bezpieczeństwa. Z analizą optymalizacyjną nieodłącznie 
związane jest także badanie wrażliwości probabilistycznej parametrów 
wytrzymatosciowych lub użytkowych na zmianę własności materiałowych, wymiarów 
lub innych parametrów opisujących konstrukcję. Zgodnie z zapisem zawartym w normie 
PN-ISO 2394, zmiana każdego czynnika mogącego zakłócić miarę niezawodności, 
powinna być połączona ze zbadaniem wpływu tej zmiany na ogólną niezawodność 


Obliczenia takie umożliwiają wskazanie „słabych? elementów konstrukcji lub 
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czynników, które w projektowaniu mają największe znaczenie i wymagają szczególnej 
uwagi projektanta. 

Jak wspomniano wcześniej, narzędzia probabilistyczne nie są obecnie 
powszechnie stosowane w analizie problemów inżynierskich. Wynika to przede 
wszystkim z ich poziomu skomplikowania uniemożliwiającego wykonywanie analizy 
niezawodnościowej bez odpowiedniego przygotowania matematycznego. 
Zautomatyzowanie obliczeń z zakresu analizy niezawodności, optymalizacji lub 
wrażliwości w stopniu zbliżonym do od dawna stosowanych i rozwijanych programów 
wykorzystujących metodę elementów skończonych prawdopodobnie jeszcze długo nie 
będzie możliwe. Wydaje się jednak istotne wykazanie, że proste i jednocześnie skuteczne 
metody probabilistyczne można z łatwością zastosować w celu wzmocnienia i 
udokładnienia procesu projektowania. 

W pracy podjęto próbę sformułowania metod łączących opis probabilistyczny z 
analizą wrażliwości odpowiedzi mechanicznej konstrukcji na imperfekcje geometryczne 
i materiałowe. Wykorzystano metody, które można scharakteryzować jako inżynierskie, 
a więc takie, które jednocześnie nie wymagają zaawansowanego oprogramowania 1 są 


łatwe w zastosowaniu. 


1.2 Teza, cel i zakres pracy 


Teza pracy 
Tezę pracy sformułowano następująco: Możliwe jest oszacowanie wrażliwości 
stanów granicznych nieliniowych modeli konstrukcji inżynierskich na imperfekcje 
geometryczne i materiałowe przy zastosowaniu metody symulacyjnej Monte Carlo wraz 
z odpowiednio  sformułowanymi technikami  redukcyjnymi,  próbkowaniem 


ukierunkowanym, metodą powierzchni odpowiedzi oraz metodą estymacji punktowej. 


Cel pracy 
W Katedrze Mechaniki Budowli, Wydziału Inżynierii Lądowej, Politechniki 
Gdańskiej od wielu lat rozwijane są metody analizy probabilistycznej. Wykorzystywane 
były one nie tylko w problemach dotyczących konstrukcji inżynierskich ale również w 
zagadnieniach z dziedziny ochrony środowiska. Zbudowano własne programy służące 
m.in.: generacji pól losowych, zastosowania metody perturbacyjnej w analizie powłok 
czy też wykorzystaniu metod: estymacji punktowej (PEM), Monte Carlo (MC), 


powierzchni odpowiedzi (RSM) i innych. Zaproponowane algorytmy i narzędzia 
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obliczeniowe są więc uniwersalne i można je stosować do analizy praktycznie każdego 
problemu inżynierskiego. Kluczowym jest jednak właściwy wybór metody a także 
przyjęcie odpowiedniego algorytmu obliczeniowego - czynniki te zależne są od 
charakteru konkretnego zadania. 

Ogólnym celem niniejszej rozprawy jest zaadoptowanie istniejących algorytmów 
do analizy wrażliwości stanów granicznych nośności i użytkowania wybranych 
konstrukcji inżynierskich na imperfekcje geometryczne 1 materiałowe. Dodatkowo, 
istniejące procedury wzbogacono o własne narzędzia, których głównym zadaniem jest 
uproszczenie i przyspieszenie procesu obliczeń. Przedstawione w rozprawie przykłady 
ograniczono do prostych układów prętowych, ilustrujących zagadnienie. 

W pracy pod pojęciem wrażliwość konstrukcji będzie rozumiana zmiana 
parametrów nośności lub użytkowalności pod wpływem imperfekcji geometrycznych lub 
materiałowych, a także obciążenia zewnętrznego. Wychodząc z założenia, że każde 
obciążenie oraz imperfekcje mogą (a nawet powinny) być zdefiniowane jako zmienne 
losowe o określonej zmienności, zdecydowano się zbadać wpływ tej zmienności, opisanej 
w sposób probabilistyczny, na odpowiedź mechaniczną konstrukcji. Efektem końcowym 
analizy jest zatem estymacja zmiany niezawodności konstrukcji. Nie należy więc 
zaproponowanego podejścia traktować jako standardowej analizy wrażliwości, można je 


określić jako niezawodnościową lub probabilistyczną wrażliwość konstrukcji. 


Podsumowując, cele niniejszej pracy można sformułować następująco: 
« przegląd probabilistycznych metod obliczeniowych pozwalających na analizę 
niezawodnościowej wrażliwości konstrukcji, 
= dobór i adaptacja wybranych metod właściwych dla proponowanej analizy, 
= zaproponowanie własnych algorytmów optymalizujących obliczenia, 
= przetestowanie zaproponowanych algorytmów przy wykorzystaniu odpowiednio 


dobranych przykładów; 


W pracy przyjęto szereg ograniczeń i założeń: 
= praca nie obejmuje analizy wrażliwości rozumianej w standardowy sposób, 
analizowane są jedynie metody wrażliwości niezawodnościowej, 
« przykłady obliczeniowe dotyczą jedynie statyki konstrukcji, chociaż możliwe jest 


bezpośrednie rozszerzenie algorytmów o analizy dynamiczne, 
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= algorytmy sformułowano za pomocą narzędzi obliczeniowych umożliwiających ich 
bezpośrednie zastosowanie w analizie konstrukcji rzeczywistych, 
= obliczenia wykonano dla prętowych modeli teoretycznych oraz przykładów 


rzeczywistych konstrukcji inżynierskich - kratownic. 


Zakres pracy 
Pracę podzielono na 6 rozdziałów. 

W pierwszym rozdziale przedstawiono tezę pracy oraz dokonano przeglądu 
literatury. Głównym celem tego przeglądu było zdefiniowanie zakresu pracy, na tle wielu 
zagadnień związanych z zastosowaniem probabilistyki w problemach inżynierskich. 
Obszerne fragmenty poświęcono publikacji polskich autorów, a także, w węższym 
zakresie gdańskiemu ośrodkowi, który od lat rozwija tę dziedzinę mechaniki. Wydaje się 
to ważne, gdyż niniejsza praca jest twórczą kontynuacją wcześniejszych opracowań. 

W drugim rozdziale przybliżono najważniejsze pojęcia związane z 
zagadnieniami probabilistycznymi a w szczególności niezawodnością konstrukcji. 
Rozdział stanowi pewne rozszerzenie przeglądu literatury i ma zdefiniować pojęcia 
związane z treścią pracy. Między innymi omówiono poziomy metody probabilistycznych 
oraz modele losowe zmiennych projektowych. Ponadto uszczegółowiono pojęcie funkcji 
stanu granicznego. 

Rozdział trzeci omawia szczegółowo bezpośrednio wykorzystywane w pracy 
metody: Monte Carlo (MC), powierzchni odpowiedzi (RSM), estymacji punktowej 
(PEM) oraz techniki redukcyjne próbkowania warstwowego (SS) i hipersześcianu 
łacińskiego (LHS). Rozdział zamyka szerszy opis techniki próbkowania 
ukierunkowanego (TRS), która została wykorzystana w budowie autorskiego programu. 

Pojęcia i metody analizy związane z badaniem wrażliwości mechanicznej 
odpowiedzi konstrukcji na zmianę jej parametrów projektowych przedstawiono w 
rozdziale czwartym. Na podstawie literatury dokonano klasyfikacji tych metod, 
przybliżając różnorodne definicje wrażliwości oraz szeroki wachlarz ich zastosowań w 
analizie zagadnień inżynierskich. Kilka pojęć opisano w szerszy sposób, gdyż zostały 
wykorzystane w obliczeniach w dalszej części pracy. Należą do nich m.in. projektowanie 
przesiewowe oraz metoda Monte Carlo (MC) ukierunkowana na estymacje indeksów 


wrażliwości Sobola. 
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W najobszerniejszy rozdziale piątym przedstawiono kilka przykładów 
wykorzystania omówionych wcześniej metod. W każdym z modeli starano się skupić 
uwagę na wybranych zagadnieniach estymacji niezawodności i wrażliwości konstrukcji, 
aby w miarę możliwości przetestować jak największa liczbę wariantów obliczeniowych. 
Należy podkreślić, że tylko praktyka obliczeniowa może wykazać przydatność 
poszczególnych metod w obliczeniach inżynierskich. 

W pierwszym przykładzie (Rozdz. 5.1) wykorzystano standardową i powszechnie 
znaną kratownicę von Misesa. Dla takiego, obliczeniowo łatwego modelu można było 
uzyskać dokładne histogramy nośności granicznej, oddzielnie dla każdej zmiennej 
projektowej. Takie podejście nie jest możliwe w przypadku analizy rzeczywistych 
konstrukcji. Obliczono globalne oraz lokalne współczynników opisujące wrażliwość 
konstrukcji, a także wskaźniki wrażliwości Sobola. 

Kolumna Ziglera, drugi z teoretycznych modeli przedstawionych w pracy 
(Rozdz. 5.2), posłużył do wykonania analizy chmur rozproszenia wyników obliczeń i w 
efekcie scharakteryzowania wrażliwości układu. Przeprowadzono także wszechstronną 
analizę modelu przy wykorzystaniu metody Monte Carlo (MC) wersji bezpośredniej oraz 
stosując metody redukcyjne próbkowania warstwowego (SS) oraz hipersześcianu 
łacińskiego (LHS). Uzyskane wyniki porównano z obliczeniami wykonanymi za pomocą 
metody powierzchni odpowiedzi (RSM), przy wykorzystaniu tych samych, wcześniej 
wygenerowanych próbek. Podrozdział zamyka analiza przeprowadzona autorską metodą 
próbkowania ukierunkowanego (TRS). 

Najobszerniejszym przykładem przedstawionym w pracy jest analiza 
rzeczywistej wieży telekomunikacyjnej (Rozdz. 5.3). Aby wykonać niezwykle 
kompleksowe i pracochłonne obliczenia porównawcze, przyjęto jedynie dwie zmienne 
opisujące sztywność podparcia wieży. Oprócz klasycznej metody Monte Carlo (MC) 
opracowano cztery warianty obliczeń za pomocą metody powierzchni odpowiedzi 
(RSM), wykonane każdorazowo dla populacji 200 próbek. Tak jak w poprzednich 
przykładach zastosowano także metodę próbkowania ukierunkowanego (TRS). 
Przeprowadzona analiza pozwoliła na określenie wrażliwości konstrukcji z uwagi na 
zmianę parametrów jej podparcia. Obszerne wyniki wykorzystano w aproksymacji 
funkcji opisujących wrażliwość, co stanowiło podstawę do wyznaczenia odpowiednich 
liczbowych parametrów. Dodatkowo przyjęto autorski opis współczynników wrażliwości 


bliski wynikom analizy wykonanej za pomocą komercyjnego programu. Końcowym 
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efektem obliczeń jest pewna forma powierzchni wrażliwości, którą można uzyskać 
jedynie w wyniku czasochłonnych obliczeń w otoczeniu punktu projektowego. Jej 
znaczenie jest więc ograniczone. 

W ostatnim przykładzie (Rozdz. 5.4) zaprezentowano kompleksowe obliczenia 
standardowej kopuły prętowej [Sorn i inni, 2016]. W tym przypadku zastosowano 
przede wszystkim metodę estymacji punktowej (PEM). Za pomocą wstępnej analizy 
parametrycznej uzasadniono wybór tej metody. Wynikiem końcowym analizy 
wrażliwości było określenie minimalnego pola przekroju poprzecznego prętów przy 
założeniu 50-cio letniego normowego poziomu niezawodności. 

Pracę zamykają wnioski (Rozdz. 6) oraz spis literatury obejmujący 204 pozycje. 
Ponieważ literatura z zakresu zastosowań probabilistyki jest niezwykle obszerna 
przedstawiony spis ma bardzo ograniczony charakter. Skupiono się na obszarze ściśle 
związanym z prezentowanymi tematami. 

W opisie metod probabilistycznych zdecydowano sie na równoległe 
wykorzystanie nazewnictwa angielskiego. Wynika to z braku odpowiedniej terminologii 
w języku polskim lub powszechnej zgody dotyczącej tego nazewnictwa. Formułowanie 
własnych określeń wprowadziłoby wiele nieporozumień i niejednoznaczności. Dotyczy 


to także przyjętych skrótów. 


Za najważniejsze elementy pracy można uznać: 
« sformułowanie algorytmów obliczeniowych umożliwiających badanie 
niezawodnościowej wrażliwości konstrukcji, 
= zbudowanie własnego programu wykorzystującego technikę próbkowania 
ukierunkowanego, pozwalającego na właściwy dobór próbek, a w efekcie na 
znaczne przyspieszenie wykonywanych obliczeń, 
=" adaptacja wybranych metod (Monte Carlo, PEM, metoda powierzchni 
odpowiedzi), w taki sposób aby można było przeprowadzić analizę wrażliwości 
niezawodnościowej. 
Stworzone procedury obliczeniowe mają szersze, uniwersalne znaczenie, gdyż rozwijają 


metody analizy probabilistycznej. 
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1.3 Przegląd literatury 


Niezawodność konstrukcji inżynierskich jest dziedziną na pograniczu nauk 
podstawowych i stosowanych. Analiza probabilistyczna, wraz z matematycznymi 
podstawami, służy do wypracowania procedur o użytkowym charakterze, 
umożliwiających usprawnienie procesu projektowania w kierunku: obliczeniowej 
efektywności, właściwego modelowania oddziaływań, materiałów i struktur oraz 
ekonomii. Podstawy niezawodności konstrukcji jako nauki zawarte są w szeregu 
opracowań monograficznych [Hart, 1982; Thoft-Christensen i inni, 1982; Augusti i inni, 
1984; Madsen i inni, 1985; Ditlevsen i inni, 1996; Melchers, 1999, Ang i inni, 2007; 
Raizer, 2004]. Istnieje bogata literatura poświęcona analizie niezawodności wybranych 
typów konstrukcji inżynierskich, ze względu na: materiał konstrukcyjny, własności 
fizyczne i przyjęty model losowy. Zagadnienia losowej analizy konstrukcji stalowych 
opisane są w np. pracach [Biegus, 1999; Guedes-Soares i Garbatov, 2002], analizie 
niezawodności konstrukcji drewnianych poświęcone są m.in. prace [Koechler, 2007; 
Pieniak i inni, 2011], a tematyka niezawodności konstrukcji betonowych ujęta jest w 
[Nowak i Szerszen, 2002]. Wyróżnione są także przypadki niszczącego działania 
zewnętrznego: wiatru [Simiu i Scanlan, 1996], oddziaływań sejsmicznych [Bołotin, 
1961], łącznie działającego falowania morskiego i wiatru w przypadku budowli typu 
offshore, gdzie konieczna jest analiza w dziedzinie czasu [Reliability of Offshore, 1992]. 

W Polsce przez lata wiodącym na tym polu był ośrodek krakowski pod 
kierownictwem prof. Janusza Murzewskiego [Murzewski, 1970, 1989]. Wartościowe 
prace w zakresie niezawodności konstrukcji metalowych wniósł prof. Zbigniew Mendera 
[Mendera, 1987]. Pod egidą środowiska krakowskiego w latach 1967 — 1992 odbyły się 
wydarzenia naukowe prezentujące metody probabilistyczne w mechanice z udziałem 
prelegentów polskich i zagranicznych [Losowe obciążenia, 1979; Projektowanie 
konstrukcji, 1988; Podstawy projektowania, 1998]. Istotny wkład w rozwój dziedziny 
wniósł ośrodek wrocławski [Biegus, 1999; Śniady, 2000], z kolei elementy losowej teorii 
obciążeń zawarte są w opracowaniu [Engel i Sieczkowski, 1978]. Pionierskim w 
dziedzinie losowej mechaniki był także ośrodek świętokrzyski, kierowany przez prof. 
Zbigniewa Kowala [Kowal, 2003]. Należy także wymienić prace: [Sobczyk, 1973, 1991; 
Sobczyk i inni, 1996; Wilde, 1981; Skalmierski i inni, 1982; Tylikowski, 1991; Kleiberi 
inni, 1992; Hien, 2003; Stocki, 1999, 2010; Woliński i Wróbel, 2001; Kolanek, 2007; 
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Kamiński, 2001, 2013; Gwóźdź i inni, 2011]. Zastosowanie metod probabilistycznych 
propagowali w swoich pracach także Bryja i Śniady [Bryja i inni, 1988], Mazur-Śniady i 
Śniady [Mazur-Sniady i inni, 1986], Sieniawska i Śniady [Sieniawska i inni, 1990], 
Śniady, Sieniawska, Żukowski [Śniady i inni, 1998], Sieniawska, Wysocka, Zukowski 
[Sieniawska i inni, 2000]. Kleiber, Knabel i Rojek [Kleiber i inni, 2004], Kleiber, 
Siemaszko i Stocki [Kleiber i inni, 1999], Knabel [Knabel. 2004], Stocki, Kolanek, Jendo 
i Kleiber [Stocki i inni, 2001], Kowal [Kowal, 2011], Dudzik i Radoń [Dudzik i inni, 
2015], Radoń [Radoń, 2011, 2015], Rzeszut, Folta i Garstecki [Rzeszut i inni, 2018]. 
Przeglądając wymienione pozycje literatury można zauważyć, że omawiają one bardzo 
różnorodne zagadnienia, począwszy od ogólnych podstaw teoretycznych i metodologii 
rozwiązania na analizie szczególnych konstrukcji lub ich fragmentów skończywszy. 
Najbardziej znanym i najczęściej stosowanym narzędziem obliczeniowym, 
ukierunkowanym na numeryczną analizę zjawisk losowych jest metoda symulacyjna 
Monte Carlo (MC), zwana metodą prób statystycznych [Metropolis i Ulam, 1949; 
Hammersley i Handscomb, 1964; Rubinstein, 1981]. Z metodą Monte Carlo 
nierozerwalnie związany jest Polski uczony Stanisław Ulam, który wykorzystał ją w 
projekcie Manhattan w ośrodku badań jądrowych w Los Alamos [Ulam, 1961]. 
Podstawowy algorytm symulacyjny zadania inżynierskiego polega na generacji 
zmiennych podstawowych w postaci liczb losowych, wykonaniu w symulacyjnej pętli 
deterministycznego zadania projektowego oraz opisaniu wyników w formie histogramu, 
np. estymatora rozkładu prawdopodobieństwa stanu granicznego. Ważnym zagadnieniem 
numerycznym, Ściśle związanym z metodą Monte Carlo, jest generacja liczb losowych o 
zadanym rozkładzie prawdopodobieństwa [Zieliński, 1970; Gentle, 1998]. W przypadku 
istnienia odwracalnej dystrybuanty, generacja ta jest możliwa poprzez transformację 
zbioru liczb losowych o rozkładzie równomiernym z przedziału [0, 1] [Devroye, 1986; 
Wieczorkowski i inni, 1997]. Podstawowa, bezpośrednia technika losowania prostego 
metodą Monte Carlo jest najłatwiejsza do oprogramowania, jednak przy dużych 
zadaniach pochłania znaczący czas obliczeniowy. Opracowano wiele algorytmów 
redukcji wariancji, ukierunkowujących losowanie na tzw. obszar krytyczny a więc 
estymacji niezawodności lub prawdopodobieństwa awarii. Rozwiniętymi, tak 
analitycznie jak i obliczeniowo, technikami redukcyjnymi są m.in.: próbkowanie z 
funkcją ważności (Importance Sampling), próbkowanie warstwowe (Stratified 


Sampling), próbkowanie kierunkowe, (Targeted Random Sampling) czy próbkowanie z 
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zastosowaniem tzw. hipersześcianu łacińskiego (Latin Hypercube Sampling) [Melchers, 
1999; Hurtado i inni, 1998; Shields 2014, 2015]. 

Użytecznym narzędziem numerycznym w analizie probabilistycznej konstrukcji 
inżynierskich jest metoda powierzchni odpowiedzi (Response Surface Method - RSM) 
[Box i inni, 1987; Hurtado i inni, 1998]. Porównanie RSM z metodą sieci neuronowych 
a także innymi probabilistycznymi metodami można znaleźć np. w [Gomes i inni, 2004]. 
Ponawiane są próby zastosowania w zagadnieniach inżynierskich teorii zbiorów 
rozmytych (Interval Analysis) [Niczyj, 2003], jednak ich wykorzystanie w badaniach 
rzeczywistych konstrukcji wydaje się ciągle odsunięte w czasie. Taką samą uwagę można 
sformułować o metodach opartych na analizie zbiorów wypukłych [Au, 2003]. Należy 
także wspomnieć o metodzie stochastycznych elementów skończonych SFEM, jako 
rozszerzeniu klasycznej metody elementów skończonych [Stefanou 2009; Sudret i inni, 
2002]. 

Na bazie różnych probabilistycznych metod wykształciły się dyscypliny 
pochodne, tworzące szeroką tematykę niezawodności. Jedną z nich jest optymalizacja 
niezawodnościowa. Jest to rozwijająca się od kilku dekad dziedzina, formułująca zadania 
teorii projektowania w kategoriach losowych - zarówno w losowej formacji integralnych 
składników funkcji celu jak 1 warunków ograniczających. Optymalizacja 
niezawodnościowa jest obecnie rozwiniętą dziedziną, zarówno w teoretycznych 
podstawach jak i aplikacjach numerycznych z dostępnymi pakietami oprogramowania 
inżynierskiego [Thoft-Christensen i inni, 1986; Der Kiureghian, 2003]. Ważna klasa 
zagadnień losowej optymalizacji niezawodnościowej jest ukierunkowana na koszt 
konstrukcji [Rackwitz, 2002, 2003; Augusti, 2003]. 

Z optymalizacją niezawodnościową związane jest pojęcie wrażliwości 
probabilistycznej stanów granicznych konstrukcji na rozrzut zmiennych podstawowych 
[Bjerager i inni, 1989; Melchers i inni, 2004]. Często stosowane jest specjalistyczne 
pojęcie „robust design”, spotykane w zakresie projektowania urządzeń technicznych czy 
elementów konstrukcyjnych. W ogólnym ujęciu podejście to kładzie nacisk na 
uwzględnienie w analizie niezawodności poszczególnych elementów składowych 
[Kececioglu, 2003; Lagaros, 2007]. 

Metodom analizy wrażliwości poświecono wiele książek i publikacji, zarówno w 
zakresie deterministycznym jak i probabilistycznym, a także łączącym te dwa podejścia. 


W przypadku metod probabilistycznych zastosowanych w analizie wrażliwości należy 


N/N MOST 


1.3 Przeglad literatury 15 


wymienić takie przeglądowe pozycje jak np. [Saltelli i inni, 2008] lub prace formułujące 
ogólne uwagi [Fang i inni, 2014]. Definicje i klasyfikacje metod związanych z analizą 
niezawodności konstrukcji można znaleźć m.in. w [Hamby, 1994; Saltelli i inni, 2000, 
2004; Saltelli, 2002; Frey i inni, 2003; Heiselberg i inni, 2009; Cukier i inni, 1973; 
Bertrand i inni, 2015; Nguyen i inni, 2015; Woods i inni, 2015; Kleiber i inni, 1991]. 

Probabilistyczna analiza wrażliwości konstrukcji, podobnie jak analiza 
niezawodności czy optymalizacja niezawodnościowa, wykonywana jest przy 
zastosowaniu szerokiego spektrum metod. Do klasyki obliczeniowej należy metoda 
Monte Carlo [Melchers i inni, 2004; Ahammed i inni, 2006; Wang i inni 2010]. 
Stosowane są także metody redukcyjne [Feng i inni, 2010]. Do popularnych technik 
zaliczyć można również metodę rozwinięcia funkcji losowej w chaos wielomianowy 
[Crestaux i inni, 2009; Gratiet 1 inni 2016] czy programowanie liniowe [Castillo i inni, 
2008; Hansson i-inni, 2006]. Poszukiwane są oryginalne rozwiązania [Wei i inni, 2012]. 
Osobny temat stanowi analiza wrażliwości układów nieliniowych [Kong, 2008]. 

Nie istnieje jedna zunifikowana miara opisująca wrażliwość. Sobol zaproponował 
jeden z najbardziej zaawansowanych sposobów prezentacji wrażliwości w postaci 
wskaźników| Sobol, 1990, 2001], których zastosowaniu poświęcone są m. in. prace 
[Homma i inni, 1996; Sudret, 2009]. Nowe propozycje wielkości opisujących wrażliwość 
można znaleźć w [Xiao i Lu, 2017]. Wrażliwości probabilistyczna stanów granicznych 
konstrukcji na zmienność losową zmiennych podstawowych przy wykorzystaniu 
histogramów można znaleźć w [Skowronek, 2006]. 

Za pomocą metod wrażliwości probabilistycznej analizowane są praktycznie 
wszystkie typy konstrukcji. Nacisk kładziony jest jednak na układy podatne na 
imperfekcje geometryczne, np. układy prętowe, analizowane w [Bojczuk; 1999; Silicki i 
inni 2010; Stocki i inni, 2001; Valdebenito i inni, 2012; Song i inni, 2009; Ma, 2011; 
Ikedaa, 2007], czy pręty cienkościenne [Kotełko i inni, 2017; Zeinoddini i inni, 2012; 
Kala, 2005; Papadopoulos i inni 2013]. Pominięto tutaj obszerną literaturę dotycząca płyt 
i powłok. Warto wspomnieć także o innych dziedzinach jak np. mechanika pękania 
[Sprung, 2007] lub analiza wrażliwości platform wiertniczych [Rozmarynowski i inni, 
2018]. 

Równolegle z rozwojem dziedziny niezawodności wykształciła się losowa teoria 
podejmowania decyzji, obecnie ilustrowana bogatą literaturą [Benjamin i inni, 1977; 


Faber, 2002]. 
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Szybki rozwój możliwości obliczeniowych komputerów umożliwił 
oprogramowanie wielu metod probabilistycznych. Powstały takie programy jak 
STRUREL [Gollwitzer i inni, 2006], Proban [Tvedt, 2006], COSSAN [Schueller i inni., 
2006] a także Numpress Explore, opracowany w IPPT w Warszawie. Zdecydowana 
większość prac wykorzystuje jednak komercyjne programy metody elementów 
skończonych (MES), budowane są jedynie algorytmy opracowujące wyniki uzyskane za 
ich pomocą. 

Na zakończenie tego krótkiego przeglądu należy podsumować pracę ośrodka 
gdańskiego, który od wielu lat zajmuje się zastosowaniami metod probabilistycznych w 
mechanice konstrukcji. Rozwój tych prac datuje się od lat 70-tych XX w.. Zespół 
badawczy przez ponad dwie dekady kierowany był przez prof. Eugeniusza Bielewicza. 
Działania środowiska gdańskiego ukierunkowane były na generację wielowymiarowych 
pól losowych do opisu geometrycznych imperfekcji konstrukcji powierzchniowych oraz 
ich oceny niezawodności: [Bielewicz i inni, 1985, 1987, 1994, 1997, 2002; Walukiewicz 
i inni, 1995; Skowronek, 1987; Górski 2006; Orlik 1976, Wilde 1981, Knabe i inni, 1998; 
Górski i inni 2008, 2015]. Ośrodek gdański był w latach 1980 i 1985 inicjatorem i 
organizatorem konferencji “Problemy losowe w mechanice konstrukcji”, stanowiącej 
ogólnokrajowe forum poświęcone metodom losowym w mechanice [Problemy losowe, 
1980, 1985]. Wiele prac było także związanych z modelowaniem pól losowych 
opisujących różne materiały i zagadnienia, np. pękanie [Anders, 2009], ochronę 
środowiska [Jankowski i inni, 1997], mechanikę gruntów [Tejchman i inni, 2009(a), 
2011; Przewłócki i inni 2001, 2016], elementy betonowe [Tejchman i inni, 2009(b); 
Korol 2012] a także betony asfaltowe [Szydłowski i inni, 2018]. W ostatnich latach 
główny nacisk położono na zastosowaniu metody powierzchni odpowiedzi (RSM) w 
zagadnieniach inżynierskich [Winkelmann, 2013; Winkelmann i inni, 2014]. 

Niniejsza praca jest kolejnym wkładem ośrodka gdańskiego w rozwój metod 
probabilistycznych w mechanice konstrukcji. Dotyczy rozwinięcia wcześniej 
opracowanych technik w kierunku badania wrażliwości mechanicznej odpowiedzi na 


losową zmianę parametrów opisujących geometrię i definiujących materiały. 
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Rozdział 2. Niezawodność konstrukcji inżynierskich - 


podstawowe pojęcia i definicje 


2.1 Wprowadzenie 


Losowa natura zjawisk zachodzących w przyrodzie od lat stanowi przedmiot 
zainteresowania zarówno filozofii jak i nauk ścisłych. Potrzeba odczuwania ładu i 
przewidywalności w życiu codziennym stała się dla uczonych motywacją do prowadzenia 
badań nad istotą losowości oraz próbą rozpoznawania tkwiących w zjawiskach losowych 
regularności dających przełożyć się na język matematyczny. Przypadający na lata 
trzydzieste ubiegłego wieku dynamiczny rozwój dwóch dziedzin matematyki: teorii 
prawdopodobieństwa oraz statystyki matematycznej umożliwił wypracowanie metod 
badania zjawisk losowych oraz zdefiniowanie ich w języku matematycznej teorii miary. 
Wypracowane techniki wydobywania i analizowania informacji zawartych w danych 
liczbowych pochodzących z obserwacji i badań doświadczalnych stanowią cenne 
narzędzie umożliwiające stworzenie modelu matematycznego jak najwierniej 
opisującego zmienność danego zjawiska. Obecnie teoria prawdopodobieństwa jest 
rozległą i dynamicznie rozwijającą się dziedziną matematyki ściśle związaną z 
różnorodnymi zastosowaniami. [Sobczyk 2010] W odniesieniu do zagadnień mechaniki 
konstrukcji procedury obliczeniowe oparte na rachunku prawdopodobieństwa oraz 
statystyce matematycznej stanowią naukową podstawę opisu zagadnień związanych z 
niezawodnością i bezpieczeństwem konstrukcji. 

Niezależnie od ciągłego wzrostu możliwości obliczeniowych oraz stopnia 
złożoności współczesnych konstrukcji inżynierskich, główne założenia dotyczące ich 
projektowania nie uległy zmianie. Zgodnie z wytycznymi zawartymi w normie PN-EN 
1990 konstrukcję należy zaprojektować tak, aby jej nośność, użytkowalność oraz 
trwałość były należyte. Sprawdzenie czy konstrukcja spełnia powyższe wymagania 
odbywa się poprzez weryfikację tzw. stanów granicznych. Stany graniczne rozdzielają 
dopuszczalne i niedopuszczalne stany konstrukcji - ich przekroczenie może doprowadzić 
do zniszczenia układu (stan graniczny nośności) bądź utrudniać jego prawidłowe 
funkcjonowanie (stan graniczny użytkowania). 

O tym w jakim stanie znajduje się konstrukcja decyduje zbiór wielkości tzw. 


zmiennych podstawowych, do których wg PN EN 1990 zalicza się: oddziaływania, 
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wpływy środowiskowe, właściwości materiałów i wyrobów oraz dane geometryczne. 
Wszystkie z wymienionych wielkości podlegają wpływom czynników losowych. 
Podstawowym problemem, przed którym stoi projektant jest uwzględnienie losowej 
natury zmiennych w taki sposób, aby obliczenia wykonywane z ich udziałem 
odzwierciedlały możliwie jak najlepiej rzeczywisty stan konstrukcji. 

Ze względu na zmienny charakter parametrów projektowych, analiza stanu 
naprężenia i odkształcenia w konstrukcji musi się wiązać z właściwym osądem czy 
konstrukcja jest w stanie służyć wystarczająco niezawodnie, przy czym według PN ISO 
239 przez niezawodność rozumie się zdolność układu technicznego do wypełniania 
funkcji przewidzianych w projekcie, w danych warunkach eksploatacji. W zależności od 
tego w jaki sposób szacowana jest niezawodność (co przyjmiemy za jej miarę a także w 
jaki uwzględniono losową naturę zmiennych projektowych) w literaturze można znaleźć 


definicję trzech poziomów metod probabilistycznej analizy konstrukcji [Madsen, 1979]. 
2.2 Poziomy metod probabilistycznej analizy konstrukcji 
2.2.1 Metody poziomu I 


Podstawy metod poziomu I podane zostały w latach 40-tych i 50-tych XX w. 
[Freudenthal, 1947, 1956]. Ich istotą jest zapewnienie niezawodności konstrukcji bądź jej 
elementu poprzez zastosowanie tzw. częściowych współczynników bezpieczeństwa, 
które modyfikują wartości nominalne bądź charakterystyczne zmiennych projektowych 
według zaleceń przedstawionych w normach. Zgodnie z klasyfikacją przyjętą według 
normy PN-ISO 2394 zmienne podstawowe można podzielić na wyrażające efekty 
oddziaływań S oraz wyrażające nośność R. Symboliczny warunek bezpieczeństwa 
przyjęty według metody częściowych współczynników w odniesieniu do stanu granicznej 
nośności ma następującą postać: 


g(S,,R,)20 (2.1) 


przy czym S, i R, są wartościami obliczeniowymi odpowiednio: efektu działań na 
konstrukcję oraz jej odporności. Wielkości te stanowią iloczyn bądź iloraz wielkości 
nominalnych/charakterystycznych zmiennych projektowych oraz częściowych 
współczynników bezpieczeństwa (współczynniki obciążeń, współczynniki materiałowe 


współczynniki konwersji itp.). 
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W przypadku stanu granicznej użytkowalności równanie ma postać: 
Sı SC (2.2) 

gdzie C jest założonym ograniczeniem wynikającym z użytkowalności elementu. 

Metoda stanów granicznych w połączeniu z metodą częściowych 
współczynników bezpieczeństwa stanowi obecnie podstawę większości norm 
dotyczących projektowania konstrukcji. Proces obliczeń zgodnie z tą metodą przebiega 
szybko oraz nie wymaga zastosowania skomplikowanych procedur zawierających 
elementy rachunku prawdopodobieństwa. Do wad metody zaliczyć można brak 
informacji dotyczących wpływu zmiany poszczególnych współczynników 
bezpieczeństwa na niezawodność konstrukcji (metody probabilistyczne umożliwiają 
ilościową ocenę prawdopodobieństwa osiągnięcia stanów granicznych w założonym 
okresie eksploatacji) oraz fakt, iż dobór właściwych współczynników bezpieczeństwa 
wymaga przyjęcia określonych założeń i nie jest wolne od decyzji arbitralnych. [Lewicki, 
1995] 

Procedury projektowania metodami poziomu I mają charakter deterministyczny, 
jednak kalibracja częściowych współczynników przeprowadzana jest na podstawie zasad 
rachunku prawdopodobieństwa i statystyki matematycznej, dlatego metody poziomu I 


określa się często mianem półprobabilistycznych. [Skowronek, 2006] 
2.2.2 Metody poziomu II 


W szacowaniu niezawodności metodami probabilistycznymi uwzględnienie 
zmiennej natury parametrów projektowych (np. wymiarów elementów, obciążenia, 
wytrzymałości materiału itp.) uzyskuje się poprzez przedstawienie ich jako zmiennych 


losowych 4X, |. W przypadku metod poziomu II, zwanych bezrozkładowymi, nie jest 


wymagana pełna identyfikacja typu rozkładu zmiennych projektowych — wystarczającą 
podstawę obliczeń stanowią dwa pierwsze momenty tych rozkładów będące miarą 
wartości centralnej i rozrzutu. 

W procedurach metod probabilistycznych poziomu II jako miarę niezawodności 
przyjmuje się tzw. wskaźnik niezawodności 5 (określany często mianem wskaźnika 
niezawodności „drugiego momentu”). Wraz z rozwojem przybliżonych metod analizy 
niezawodności, ewolucji ulegał również wskaźnik f, przyjmując w swej nazwie 


rozszerzenie oznaczające sposób jego szacowania [Knabel, 2004]. 
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Gdy wszystkie zmienne podstawowe opisane są rozkładem Gaussa możliwe jest 


algebraiczne sformułowanie kryterium stanu granicznego, czyli tzw. funkcji stanu 
granicznego g (X), której argumentami są zmienne podstawowe. W przypadku liniowej 


funkcji stanu granicznego jest ona automatycznie zmienną o rozkładzie normalnym, 
natomiast estymacja prawdopodobieństwa awarii jest standardowym matematycznym 
zadaniem w zakresie zmiennych Gaussa. Nieliniową funkcję stanu granicznego można 
linearyzować np. w punkcie wartości średnich, dochodząc do podobnego rezultatu. 
Powyższe działania prowadzą do tzw. wskaźnika niezawodności Cornella [Cornell, 1969] 


— ilorazu estymatorów wartości średniej i odchylenia standardowego funkcji stanu 


granicznego g(X). Ogólne sformułowanie wskaźnika Cornella obciążone jest wadą — 


przy różnych matematycznych formach funkcji g(X) opisujących to samo zadanie 


można otrzymać różne wskaźniki niezawodności, a tym samym różne 
prawdopodobieństwa awarii. Rozwiązanie nie jest więc niezmiennicze [Putresza i inni, 
1995]. 

Propozycję rozwiązania niezmienniczego względem transformacji funkcji stanu 
granicznego podali Hasofer i Lind [Hasofer i inni, 1974]. W koncepcji tej nieliniowa 
funkcja stanu granicznego linearyzowana jest w tzw. punkcie projektowym, o 
najmniejszej odległości od początku układu w przestrzeni zmiennych podstawowych 
standaryzowanych. Kolejne koncepcje wskaźników niezawodności podali Rackwitz i 
Fiessler [Rackwitz i inni, 1978] oraz Ditlevsen [Ditlevsen, 1979]. 

W zakresie metod poziomu II wyróżniamy dwa podejścia szacowania 
wskaźników niezawodności: FORM (First-Order Reliability Methods), w którym funkcja 
stanu granicznego aproksymowana jest funkcją liniową oraz SORM (Second-Order 


Reliability Methods), o aproksymacji kwadratowej [Melchers, 1999; Nowak, 2000] 
2.2.3 Metody poziomu III 


Analiza niezawodności przeprowadzona metodami poziomu III pozbawiona jest 
wszelkich uproszczeń, przez co, w przypadku rzeczywistych problemów inżynierii, nie 
jest możliwe zastosowanie zamkniętych rozwiązań analitycznych. Dotyczy to zarówno 
zmiennych projektowych (w metodach poziomu III wymagana jest pełna wiedza 


odnośnie ich rozkładu) jak również samej miary niezawodności, którą w tym przypadku 
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stanowi wprost prawdopodobienstwo, ze w założonym procesie eksploatacji nie zostaną 
osiągnięte stany graniczne. 

Do zadań metod poziomu III, poza wymiarowaniem i weryfikacją niezawodności 
konstrukcji a także kalibrowaniem miar niezawodności metod niższych poziomów należy 
również optymalizacja zadań niezawodności ustrojów. 

Należy podkreślić, że w literaturze można znaleźć inne definicje poziomu III 


niezawodności, np. [Murzewski, 1989]. 
2.3 Podstawowe pojęcia i definicje 


Celem  probabilistycznej analizy konstrukcji jest oszacowanie 
prawdopodobieństwa awarii (rozumianej jako przekroczenie stanu granicznego) na 
podstawie przyjętego matematycznego modelu, w którym wybrane bądż też wszystkie 


zmienne projektowe przyjęte są jako zmienne losowe . Podstawowe parametry losowe 
grupowane są w wektorze kolumnowym X ={X,,X),...,X,} zwanym wektorem 
losowym podstawowych zmiennych projektowych. Wektor X przyjmuje wartości w n- 
wymiarowej przestrzeni liczb rzeczywistych R”. Pojedynczy zbiór wartości 
poszczególnych zmiennych podstawowych nazywamy realizacją zmiennej losowej X i 
oznaczamy jako x ={x,,X,,....x,}. 

Stany graniczne w analizie niezawodności wyrażone są poprzez tzw. uogólnioną 
funkcję stanu granicznego g (X), która dzieli przestrzeń zmiennych następująco: 

g(X)<0: obszar awarii Q, 


g(X)=0: powierzchnia graniczna (2.3) 
g(X)>0: obszar bezpieczny Q, 


Ilustrację powierzchni granicznej oraz obszarów: bezpiecznego i awarii dla przypadku 


dwuwymiarowego przedstawiono na rys 2.1. 


Na podstawie zapisanego warunku stanu granicznego prawdopodobieństwo awarii 


wyraża się następująco: 


P;=Pp(Q,)=p(s(X)<0) (2.4) 
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AX, 


LJ OBSZAR AWARII 
g(X)<0 


POWIERZCHNIA GRANICZNA 
g(X)=0 


OBSZAR BEZPIECZNY 
g(X)>0 


Rysunek 2.1 Ilustracja powierzchni granicznej oraz obszarów: bezpiecznego (Q, ) i awarii (Q, ) - 


przypadek dwuwymiarowy 


W przypadku gdy dla danego rozkładu prawdopodobieństwa istnieje funkcja gęstości to 


prawdopodobieństwo awarii można zapisać w postaci całki: 


p= | f(x)dx (2.5) 
Qp 


Wyznaczenie całki wyrażonej wzorem (2.5) może być kłopotliwe ze względu na to, że: 


« funkcja gęstości prawdopodobieństwa f(x) może nie być dobrze znana (dane 


statystyczne lub eksperymentalne mogą być niekompletne), 


= funkcja graniczna g(X ) może zawierać w sobie pewne niedokładności związane 


z przyjętym modelem zachowania się konstrukcji, 

=" całkowanie numeryczne wyrażenia 2.5 w przypadku dużej liczby zmiennych 
losowych (n>5) jest bardzo trudne, a czasem wręcz niemożliwe do 
przeprowadzenia. 

Ze względu na wyżej wymienione ograniczenia szacowania prawdopodobieństwa 
awarii p, szerokie zastosowanie znalazła metoda symulacyjna Monte Carlo, która 
umożliwia wykonanie całkowania wielowymiarowego poprzez zastosowanie 
odpowiednio dużej liczby symulacji polegających na generowaniu, zgodnie z przyjętym 
rozkładem prawdopodobieństwa, zmiennych projektowych o charakterze losowym. 
Metoda Monte Carlo została szczerzej opisana w dalszej części pracy. 


Alternatywą dla prawdopodobieństwa awarii p, , przy określaniu bezpieczeństwa 


konstrukcji, jest wskaźnik niezawodności Ø . Poniżej przedstawiono wybrane definicje 


miar niezawodności wykorzystane w niniejszej rozprawie. Szczegółowy opis tych miar 


można znaleźć m.in. w [Winkelmann, 2013]. 
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2.4 Wybrane miary niezawodności 
Przyjmując założenie, że całkowita wytrzymałość materiałów konstrukcji R oraz 
suma obciążeń S na nią działających są zmiennymi losowymi o rozkładzie normalnym, 


prawdopodobieństwo awarii p, można opisać następująco: 


- P(R-S<0)= P(u <o)a| 2 (2.6) 
ŚW o m] | 
gdzie: M jest zapasem bezpieczeństwa (także opisanym rozkładem normalnym), ©( -) 


jest funkcją dystrybuanty rozkładu normalnego o zerowej wartości średniej i 
jednostkowym odchyleniu standardowym, E[M] reprezentuje wartość oczekiwaną 
(średnią) zapasu bezpieczeństwa, a D[M ] — jego odchylenie standardowe. 

Wskaźnik niezawodności Ø, określany jest jako odwrotność współczynnika 


zmienności zapasu bezpieczeństwa [Cornell, 1969]: 


| _ DM) 
p= EM] (2.7) 


gdzie: v,, to współczynnik zmienności zapasu bezpieczeństwa. 


Interpretację graficzną wskaźnika niezawodności Cornella przedstawia rys. 2.2. 


flea) 
M>0 


stan | stan 
awarii | bezpieczny 


| Oy | Gu | M 


Box, Bou 


Rysunek 2.2 Graficzna interpretacja wskaźnika niezawodności wg Cornella 
[Winkelmann, 2013] 


Związek pomiędzy wskaźnikiem niezawodności Ø a prawdopodobieństwem zniszczenia 


można więc przedstawić następująco: 


B=" (p,) (2.8) 
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gdzie o'(-) jest funkcją odwrotną dystrybuanty rozkładu normalnego o zerowej 


wartości średniej i jednostkowym odchyleniu standardowym. 

Jeżeli zmienne R i S są nieskorelowane, to wskaźnik niezawodności /8 można 
zapisać w postaci: 

pz maa (2.9) 
|Op +05 

gdzie Ll, to wartość oczekiwana zmiennej losowej opisującej wytrzymałości materiałów, 
Hs to wartość oczekiwana zmiennej losowej obciążenia, o, oraz o, to odchylenia 
standardowe tych zmiennych. 


W szczególnych przypadkach funkcję stanu granicznego g(x) można 


zdefiniować jako n-wymiarową płaszczyznę: 
g(x)=a +) a,x; = a, ta'x (2.10) 
i=l 


gdzie: a, jest wyrazem wolnym, n wymiarem płaszczyzny, a= {a,,q,, ... a,; wektorem 
kolumnowym współczynników kierunkowych n-wymiarowej płaszczyzna, a 


X= {X),X>, ... X,} wektorem zmiennych losowych. 

W takim przypadku wskaźnik niezawodności będzie przedstawiał wzór: 
_ at a’ E[x] 

-Aa"Ca 


gdzie: E[x] jest wektorem wartości oczekiwanych zmiennej losowej x,a C, jest 


B (2.11) 


macierzą kowariancyjną zmiennej losowej x. 


Interpretację graficzna indeksu niezawodności przedstawiono na rys. 2.3 oraz rys. 2.4. 


g(x)<0 
Q, — obszar awarii 


g(x)—0 
obszar graniczny 


Q; 


obszar bezpieczny 


Rysunek 2.3 Graficzna interpretacja wskaźnika niezawodności Cornella - liniowa n-wymiarowa 
płaszczyzna stanu granicznego [Winkelmann, 2013] 
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g(x)<0 
Q, — obszar awarii 


g(x)=0 
Q, — obszar graniczny | 
aan 


Q, — obszar bezpieczny 


Rysunek 2.4 Graficzna interpretacja wskaźnika niezawodności — nieliniowa n-wymiarowa 
powierzchnia stanu granicznego [Winkelmann, 2013] 


Jedną z możliwych interpretacji wskaźnika niezawodności, zaproponowaną przez 


Hasofea i Linda w 1974 roku [Hasofer i inni, 1974] jest w tym przypadku: 


By, = min 8(x)= min J(x-E[x])” C(x- Ex], xe Q (2.12) 
g(x)=0 g(x)=0 
Punkt x leżący na powierzchni granicznej nazywa się punktem projektowym 
(obliczeniowym). 
Jeżeli powierzchnia graniczna zdefiniowana jest jako płaszczyzna n — wymiarowa 
(rys. 2.4), to wartości liczbowe wskaźnika niezawodności Hasofera — Linda ,, oraz 
wskaźnika niezawodności Cornella Ø są jednakowe. Wskaźnik niezawodności Hasofera 
— Linda jest więc uogólnieniem wskaźnika niezawodności Cornella dla nieliniowych 
powierzchni granicznych. 
Poszukiwanie wskaźnika niezawodności Hasofera — Linda /8,, można 
zdefiniować jako proces optymalizacji z jedną funkcją celu. Celem optymalizacji jest 
oczywiście znalezienie właściwego punktu projektowego x. Rozwiązania tego 


problemu dostarcza wiele algorytmów iteracyjnych, jednak nie można mieć pewności czy 


algorytmy te będą zbieżne we wszystkich przypadkach funkcji stanów granicznych. 
2.5 Losowe modele podstawowych zmiennych projektowych 


Przyjęcie właściwego modelu losowego 1 związanych z nim zmiennych 
podstawowych w procesie analizy i projektowania zależy od dostępności takich 
czynników jak: bazy rzeczywistych danych pomiarowych, programów obliczeniowych 


oraz podstaw teoretycznych dostępnych badaczowi na gruncie teorii 
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prawdopodobieństwa. Proces budowania modelu losowego w ogromnej mierze zależy 
również od celu działania, którym może być praca badawcza o charakterze naukowym 
bądź też zastosowanie inżynierskie. W przypadku analiz inżynierskich dąży się do 
maksymalnego uproszczenia stosowanych algorytmów. Podstawowymi kategoriami w 
modelowaniu losowym parametrów są: 
= zmienna losowa jednowymiarowa - wielkość jednowartościowa (pojedynczy 
wynik pomiaru), funkcja zdarzenia elementarnego; 
= zmienna losowa wielowymiarowa (wektor losowy) - kilka wyników pojedynczego 
pomiaru (np. granica plastyczności, wytrzymałości, moduł sprężystości w 
statycznej próbie rozciągania); 
= jednowartościowy proces stochastyczny (losowy) - pojedynczy, losowy wynik 
doświadczenia, w funkcji czasu - funkcja zdarzenia elementarnego i czasu [Augusti 
i inni, 1984; Ditlevsen i inni, 1996; Bucher, 2009]; 
= wielowymiarowy proces stochastyczny (losowy) - wektor składający się z 
procesów składowych, funkcji zdarzenia elementarnego i czasu [Wen i inni, 1981]; 
= pole losowe - zmienna losowa zdefiniowana w każdym punkcie dwu- lub 
trójwymiarowej przestrzeni. Przykładem są losowe odchyłki płaszcza zbiornika 
walcowego [Adler, 1981; Vanmarcke, 1983], 
= czasoprzestrzenne pole losowe — funkcja zdarzenia elementarnego i czasu, 


zdefiniowana w każdym punkcie dwu- lub trójwymiarowej przestrzeni. 


Warto podkreślić, że modele procesów losowych stosowane są w zagadnieniach 
związanych z dominującym wpływem funkcji czasu (zagadnienia dynamiki 
stochastycznej związanej z losowymi drganiami) na mechaniczną odpowiedź układu 
[Augusti, 1984; Bucher, 2009]. Z kolei modele pól losowych, o rozkładzie przestrzennym 
stosuje się w przypadku znaczącego wpływu geometrii układu [Vanmarcke, 1983; Adler, 
1981; Górski i inni, 2011]. 

Pierwszym etapem analizy niezawodności konstrukcji jest zdefiniowanie 
zmiennych podstawowych problemu. Na podstawie statystycznej bazy danych, w 
zależności od postawionego zadania, można poszukiwać estymatorów wartości średniej, 
wariancji, momentów wyższych rzędów, czy też estymatora funkcji gęstości 
prawdopodobieństwa. Ten ostatni zawiera najwięcej informacji, gdyż umożliwia 


przypisanie danej zmiennej podstawowej znanej funkcji gęstości prawdopodobieństwa, a 
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tym samym odpowiada poziomowi III analizy niezawodnosci. Analiza na poziomie 
estymatorów wartości średniej i wariancji prowadzi do poziomu II, wykorzystującego 
wskaźniki niezawodności. 

Odpowiednie przypisanie właściwego typu rozkładu prawdopodobieństwa 
poszczególnym zmiennym podstawowym w problemach inżynierii było i jest celem 
wielu prac badawczych. Wykorzystywane są w tym celu zazwyczaj bazy danych 
pomiarowych lub analiza statystyczna wyników badań eksperymentalnych. Próba 
podsumowania tych prac zawarta jest między innymi w książkach Nowaka i Collinsa 


[Nowak i inni, 2000] oraz Melchersa [Melchers, 1999]. Przykładowo: 


= w modelowaniu obciążeń stałych oraz długotrwałych składników obciążeń 
zmiennych adekwatną jest zmienna o rozkładzie normalnym [Nowak i inni, 2000]; 

=" obciążenia zmienne modelowane są zazwyczaj rozkładem gamma [Corotis i inni, 
1977], przy czym obciążenie maksymalne można modelować rozkładem 
ekstremów typu I (Gumbela) [Gumbel, 1958; Chalk i inni, 1980; Kotz i inni, 2000]. 
Odrębne podejście, ściśle powiązane z eksperymentami w przypadku zmiennego 
obciążenia mostów, zaproponowali Nowak i Hong [Nowak i inni, 1991]; 

=" obciążenie wiatrem - w uproszczonym modelu zmiennej losowej prędkości wiatru 
przypisywany jest rozkład ekstremów typu I (Gumbela) [Ellingwood i inni, 1999]. 
W pełnej analizie obciążenie wiatrem winno być modelowane procesem 
stochastycznym lub czasoprzestrzennym polem losowym [Żurański, 1978; Simiu i 
inni, 1996; Nowak i inni, 2000]; 

=" obciążenie śniegiem modelowane jest rozkładem logarytmiczno-normalnym, 
rozkładami ekstremów typu I (Gumbel) lub typu II (Frechet) [Ellingwood i inni, 
1983]. Kalibracja obliczeniowych wartości obciążenia śniegiem na podstawie 
odpowiednich statystyk przedstawiona jest w pracach [Ellingwood, 1981; 
Ellingwood i inni, 2017]; 

u obciążenie sejsmiczne (dynamiczne) a ściślej przyspieszenie podłoża gruntowego, 
jako skomplikowany proces losowy o długim okresie powrotu (rzędu 50 lat), z 
możliwym uproszczonym opisem quasistatycznym, może być modelowane 
rozkładem ekstremów typu II [Bołotin, 1961; Augusti i inni, 1984]; 

= imperfekcje geometryczne i materiałowe modelowane są na podstawie 


statystycznej bazy pomiarów odchyłek wymiarów geometrycznych; w przypadku 
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układów powierzchniowych konieczny jest model pola losowego odchyłek 
[Biegus, 1999]; 

=" obciążenie wiatrem i falowaniem morskim (obiekty offshore) modelowane są jako 
wielowymiarowe procesy losowe [Wen i inni, 1981; Simiu i Scanlan, 1996]; 

= parametry wytrzymałościowe (np. granica plastyczności stali) — modelowane są 
rozkładem logarytmiczno- normalnym lub rozkładem ekstremów typu III [Weibull, 
1951; Thoft-Christensen i inni, 1986]. Rozkład Weibulla użyteczny jest w 
modelowaniu wytrzymałości zmęczeniowej oraz w zagadnieniach mechaniki 
pękania, 

Odrębnym i trudnym w swej ilościowej analizie jest problem szacowania 


prawdopodobieństwa błędu ludzkiego [Nowak i inni, 2000] 
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Rozdział 3. Wybrane metody analizy niezawodności 


konstrukcji 


3.1 Metoda Monte Carlo MC 
3.1.1 Ogólne założenia metody 


Oficjalne narodziny metody Monte Carlo miały miejsce w 1949 roku, kiedy to w 
„Journal of the American Statistical Association” zamieszczono publikację Nicolasa 
Metropolisa i Stanisława Ulama pt. : „The Monte Carlo Method”. Rozważania zawarte w 
pracy dotyczyły wykorzystania liczb pseudolosowych do symulacji rozpraszania i 
absorpcji neutronów w ramach prowadzonego w latach czterdziestych ubiegłego stulecia 
programu budowy pierwszej bomby atomowej [Chłopek, 2009]. 

Wraz z intensywnym wzrostem mocy obliczeniowych komputerów, 
przypadającym na lata 70. XX wieku, metoda Monte Carlo zyskała na znaczeniu i 
znalazła zastosowanie w szerokiej gamie dziedzin nauk technicznych, w tym również w 
zagadnieniach dotyczących niezawodności układów konstrukcyjnych. Główną ideą 
podejścia Monte Carlo jest modelowanie matematyczne zjawisk poprzez zastosowanie 
algorytmów wykorzystujących generację liczb losowych. Podstawową zaletą metody jest 
możliwość zastosowania jej do procesów, których znaczna złożoność uniemożliwia 
przewidzenie ich wyników za pomocą podejścia analitycznego. 

Algorytm symulacji Monte Carlo, w dużym uproszczeniu można przedstawić 

następująco: 

Etap l: | wybór zmiennych projektowych opisanych losowo oraz wygenerowanie ciągu 
liczb dla każdej zmiennej zgodnie z przyjętym dla niej rozkładem 
prawdopodobieństwa, 

Etap 2: sformułowanie deterministycznej miary niezawodności układu, 

Etap 3: rozwiązanie dla każdej realizacji zmiennych podstawowych formuły logicznej 
stanowiącej miarę niezawodności układu oraz zapisanie rezultatu liczbowego, 

Etap 4: interpretacja zbioru wyników - rezultatem symulacji jest unormowany 
histogram będący estymatorem gęstości rozkładu prawdopodobieństwa 
zmiennej wynikowej. Na jego podstawie szacowane są charakterystyki 


statystyczne zmiennej wyjściowej. 
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Efektywne wykorzystanie symulacji Monte Carlo w analizie niezawodności 
wymaga po pierwsze przyjęcia modeli losowych parametrów projektowych w taki 
sposób, aby jak najwierniej opisywały daną wielkość, a po drugie opracowania 
efektywnej metody ich symulacji (wiarygodne oszacowanie prawdopodobieństwa awarii 


powinno być możliwe dla jak najmniejszej liczby symulacji). 
3.1.2 Teoretyczne podstawy metody Monte Carlo 


Bezpośrednia symulacja Monte Carlo polega na wygenerowaniu zgodnie z 


założoną dla rozpatrywanego problemu funkcją gęstości rozkładu prawdopodobieństwa 
f(x) zbioru realizacji x,, i=l..,N zmiennych losowych, a następnie 
zweryfikowaniu czy dana realizacja należy do obszaru awarii Q, czy też nie [Knabel, 


2004] . Celem dokonanego sprawdzenia jest określenie funkcji charakterystycznej zbioru 


Q , zdefiniowanej w następujący sposób: 


l gdy xeQ, 
I = ; 
sj (x) F gdy xQ, (31) 


Funkcja 7 o; (X) jest zmienną losową o rozkładzie dwupunktowym: 
pila, (X)=1=p;. a (X)=0]=1-p, (3.2) 
Wartość średnią i wariancję zmiennej 7 T (X) opisują następujące zależności: 


To, (X)=E(1q, (X))=1-p,+0-(1-p,)=p; (3.3) 


2 


Var| lo, (x)]=£| (Jo, (x) |-(E[a, (x)|) =p,(l-p,) (3.4) 


W metodzie Monte Carlo wyznaczenie stosunku liczby realizacji nalezacych do obszaru 
awarii do całkowitej liczby realizacji N umożliwia estymacje prawdopodobieństwa 


awarii p,p. Wielkość ta równa jest estymatorowi wartości średniej funkcji 


charakterystycznej zbioru Q, : 


By =o aya, (X (3.5) 


Błąd estymatora prawdopodobieństwa awarii wyraża się poprzez jego wariancję 
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n 1 LQ 1 1 
Var(By)= r 2er (la, (X:))= 5z; (-p,)=>Pr(l-P/) (3.6) 


Z powyższego zapisu wynika, że wraz ze wzrostem liczby realizacji N zmniejsza się 
błąd estymatora prawdopodobieństwa awarii p,. Miarą błędu jest tzw. współczynnik 


zmienności estymatora: 
V. = - (3.7) 


Powyższe równanie nasuwa następujące wnioski: po pierwsze współczynnik 
zmienności estymatora zmierza do zera wraz ze wzrostem liczby symulacji N, a po 
drugie nie zależy on od rozmiaru wektora zmiennej losowej X, a jedynie od wielkości 
szacowanego prawdopodobieństwa awarii. W przypadku rzeczywistych konstrukcji 
inżynierskich, dla których prawdopodobieństwo awarii oscyluje w granicach 107 + 10%, 


oznacza to, że wymagana liczba symulacji powinna mieścić się w przedziale 


N=10°+10°. W praktyce stosowany jest szereg technik, mających na celu 
zredukowanie liczby symulacji przy jednoczesnym zachowaniu ustalonego poziomu 
współczynnika zmienności wyników. Do technik tych zaliczamy między innymi: 
próbkowanie ukierunkowane (TRS - Targeted Random Sampling), próbkowanie 
warstwowe (SS - Stratified Sampling), próbkowanie hipersześcianu łacińskiego (LHS - 
Latin Hypercube Sampling), próbkowanie wagowe (IS - Importance Sampling), 
próbkowanie chybił-trafił (Hit-or-Miss Sampling), próbkowanie rosyjskiej ruletki 
(Russian Roulette Sampling) oraz wiele innych. W pracy ograniczono się do opisu trzech 


pierwszych z wymienionych technik. 


Do niewątpliwych zalet metody Monte Carlo należą: 

= możliwość zastosowania w zaawansowanych problemach, także nieliniowych (o 
dużej liczbie zmiennych lub skomplikowanej funkcji stanu granicznego); 

= możliwość pominięcia skomplikowanych wzorów i procedur stosowanych w 
innych metodach szacowania niezawodności; 

=" brak wymogu użycia skomplikowanego specjalistycznego oprogramowania, 

=" w przypadku zadań o niewielkim czasie wykonania pojedynczej symulacji metoda 
Monte Carlo w klasycznej postaci może być stosowana do weryfikacji wyników 


analizy niezawodności otrzymanych innymi metodami (np. w przypadku 
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stosowania metod FORM lub SORM metoda Monte Carlo umożliwia wykrycie 
znacznych błędów wynikających z nieuwzględnienia wielokrotnych punktów 


projektowych). 


Wadami metody są m. in.: 

u duża czasochłonność obliczeń; 

= brak pewności co do otrzymanego wyniku - dokonując nawet bardzo dużej liczby 
iteracji otrzymane wyniki zawsze będą stanowiły jedynie próbę wszystkich 
możliwości, 

= rezultaty otrzymane za pomocą metody Monte Carlo są zależne od jakości 
generatora liczb pseudolosowych (w przypadku pojedynczych zmiennych 
losowych) oraz programów symulujących pola losowe (w przypadku zagadnień 


wielowymiarowych). 
3.2 Metoda powierzchni odpowiedzi (Response Surface Method RSM) 
3.2.1 Ogólne założenia metody 


Metodą powierzchni odpowiedzi (RSM) w ogólności określa się zbiór 
matematycznych i statystycznych technik, mających na celu określenie interakcji 
pomiędzy wielkością objaśniającą interesujące nas zjawisko a zmianą wartości 
składających się na nie czynników. 

W pełnym procesie wykorzystania metody powierzchni odpowiedzi można 
wyróżnić trzy etapy. Pierwszym z nich jest wyodrębnienie spośród wszystkich badanych 
czynników tych, które w istotny sposób wpływają na powierzchnię odpowiedzi. Wstępna 
interpretacja interakcji pomiędzy badanymi czynnikami a ich wpływem na odpowiedź 
układu nazywana jest fazą zerową metody powierzchni odpowiedzi. Odpowiednie 
zaplanowanie eksperymentu daje możliwość przeglądu przestrzeni realizacji zmiennych 
zadania co umożliwia uzyskanie możliwie jak największej ilości informacji o badanym 
procesie, przy jednoczesnej minimalizacji kosztów obliczeń. Kolejne stadium metody 
skupia się na takim wymodelowaniu odwzorowania układu aby jak najlepiej przybliżało 
rzeczywiste zachowanie konstrukcji pod wpływem zmiany rozpatrywanych czynników. 
Ostatni etap metody związany jest z poszukiwaniem założeń optymalizacji mających na 


celu minimalizację liczby realizacji zmiennych zadania przy jednoczesnym zachowaniu 
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pożądanej wartości odpowiedzi konstrukcji [Bogucki 1 inni, 2003], [Montgomery, 1997], 
[Hill i inni, 1966]. 


3.2.2 Teoretyczne podstawy metody powierzchni odpowiedzi 


W analizie niezawodności zależność y=g(X) opisująca stan konstrukcji ma 
zwykle postać niejawną. Głównym celem metody powierzchni odpowiedzi jest 
zastąpienie powyższej relacji adekwatną aproksymacją y=g(x ), umożliwiającą 
efektywne oszacowanie prawdopodobieństwa awarii p,. Funkcja y(x) jest 


odzwierciedleniem relacji odpowiedzi uktadu w zaleznosci od zalozonych w zadaniu 


losowych zmiennych projektowych. 


a 


y(x)=/f(24,%,..,%,)+E (3.8) 


W powyższej formule x = |x; x3; x, | oznacza wektor realizacji zmiennych 


projektowych, natomiast € jest statystycznym błędem obliczeniowym, uwzględniającym 
fakt, że pewne źródła zmienności nie zostały zawarte w funkcji f(x). Błąd ten 
charakteryzuje się rozkładem normalnym o średniej 4, = 0 

Funkcja f(x) zwykle ma postać wielomianu pierwszego bądź drugiego stopnia. 
Jeżeli obszar, na którym poszukuje się odpowiedzi układu jest wąski, bądź charakteryzuje 
się niskim współczynnikiem skośności, wówczas powierzchnia odpowiedzi może być 
aproksymowana wielomianem stopnia pierwszego - mówimy wtedy o modelu 


pierwszego rzędu . Jego ogólna postać jest następująca: 
y(x)=Byt+ > Bx, +e (3.9) 
i=l 


Jezeli krzywizna powierzchni odpowiedzi jest na tyle znaczaca, ze nie mozna jej 
przybliżyć wielomianem stopnia pierwszego, wówczas do jej opisu zaleca się użycie 
wielomianu stopnia drugiego- model taki nazywa się modelem drugiego rzędu i zapisuje 
następująco: 
d(x) = B+ > Bx, + PB +Y YB xx +e (3.10) 
i=l i=l i<j j=2 
W celu dopasowania współczynników B wielomianu (nazywanych także 


współczynnikami regresji) stosuje się najczęściej metodę najmniejszych kwadratów. 


N/N MOST 


34 Rozdział 3 Wybrane metody analizy niezawodności konstrukcji 


Na ogół wszystkie rozwiązania metodą powierzchni odpowiedzi wykonane są przy 
użyciu jednego z powyższych modeli lub ich kombinacji. Przyjęty model powinien 
najdokładniej jak to możliwe odzwierciedlać wszelkie występujące w odpowiedzi układu 
nieliniowości (jeśli takie w ogóle istnieją). 

Mało prawdopodobnym jest aby model wielomianu był wiarygodnym 
przybliżeniem wyjściowej odpowiedzi układu na całej przestrzeni zmiennych, natomiast 
na ograniczonym obszarze pozwala on osiągnąć zadowalające dopasowanie 
[Montgomery, 2012]. 

Wykorzystanie metody powierzchni odpowiedzi w szacowaniu niezawodności 
stanowi skuteczne narzędzie w pozyskiwaniu interesujących informacji o badanym 
modelu przy jednoczesnym niskim nakładzie obliczeń numerycznych. Stosując tę metodę 
należy mieć jednak na uwadze, iż bazuje ona jedynie na aproksymacyjnie wygenerowanej 
powierzchni odpowiedzi układu a zatem wyniki uzyskane przy jej użyciu są przybliżone 
i jednocześnie zależne od poprawności przyjętej powierzchni. 

Gładka natura przyjętej aproksymacji wielomianem pozwala wyeliminować 
trudności obliczeniowe jakie pojawiają się w przypadku dyskretnego zbioru punktów 
obliczeniowych. Ma to znaczenie np. w obliczeniach wykorzystujących pochodne funkcji 
- założone aproksymowane odwzorowanie odpowiedzi układu umożliwia wyznaczenie 
pochodnych w miejscach, w których wyjściowo byłaby ona trudna do oszacowania bądź 
w ogóle nie istniała [Sobieski i inni, 1997]. 

Ograniczeniem metody jest brak możliwości stworzenia jednej odpowiedzi 
układu. Problem ten może pojawić się już na etapie planowania doświadczenia (doborze 
dyskretnych punktów na podstawie których przeprowadza się aproksymację). Liczba tych 
punktów stanowi obliczeniowy problem w przypadku skomplikowanych konstrukcji. Do 


ich ograniczenia stosuje się techniki redukcji wariancji. 
3.3 Metoda estymacji punktowej (Point Estimate Method PEM) 


Głównym założeniem metody estymacji punktowej (Point Estimate Method), po 


raz pierwszy zaproponowanej przez Rosenblueth’a w 1975 roku, jest zastąpienie ciągłej 


zmienne losowej X o funkcji gęstości rozkładu f,(x) zmienną losową dyskretną 


składającą się z N impulsów o rozkładzie prawdopodobieństwa opisanym następująco: 
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p, (x)= >) 6(x-x,)p(x,) (3.11) 


W powyższej formule 6O(x—x,) oznacza delte Diraca, natomiast p(x,) sa 
prawdopodobieństwami przypisanymi pewnym ustalonym punktom x,. Wzór 3.11 
opisuje nieskończenie dużą liczbę rozkładów. W metodzie PEM dąży się do znalezienia 
możliwie małej liczby M punków, dla której momenty rzędu k wyznaczone dla p, (x) 
będą takie same jak funkcji gęstości f,(x). Liczba punktów N wybranych do 
oszacowania zależy od rzędu k momentów probabilistycznych uwzględnionych w 
analizie. Najczęściej wybierane są dwa punkty, które w dość dokładny sposób 
umożliwiają oszacowanie wartości średniej oraz wariancji funkcji losowej. 


Porównując trzy pierwsze momenty probabilistyczne określone dla ciągłej oraz 


dyskretnej zmiennej losowej otrzymamy: 


m, = |», (3)dr=)xp(s) 


—co 


o = f em) a aan) z0 (3.12) 


gdzie m, oznacza wartość średnią, o,- odchylenie standardowe a y jest 
współczynnikiem asymetrii zmiennej losowej X. 
W przypadku gdy istnieje zmienna Y będąca pewnym przekształceniem zmiennej 


X takim, że: y = g(x), możliwe jest wyznaczenie w prosty sposób momentów rozkładu 
P,(y), przybliżających funkcję rozkładu zmiennej ciągłej f(y). 
Wartość oczekiwaną k -tej potęgi dyskretnej zmiennej Y zastępującej ciągłą zmienną 


J,(y) można przybliżyć w następujący sposób: 


e[r]-Że' (x,) p(x,) (3.13) 


W najbardziej ogólnym przypadku, kiedy funkcja Y jest funkcją zmiennej losowej X, 
której wartość średnia, odchylenie standardowe oraz współczynnik skośności są znane, 


przy zastąpieniu ciągłej zmiennej losowej X zmienną dyskretną, reprezentowaną przez 


N/N MOST 


36 Rozdział 3 Wybrane metody analizy niezawodności konstrukcji 


prawdopodobieństwa p(x ) =P. oraz p(x,)=P, (zwane wagami) przypisane punktom: 
x =X_ 1x, =x, zapisać można następujące równania: 


P +P =1 
m =xP +x P 


o: =(x.-m,)'P.+(x,—m,) P 223 


yo; =(x.—m,) P.+(x,-m,) P, 
Przekształcając powyższy układ równań otrzymujemy wartości współrzędnych punktów 
oraz przypisane im wagi: 
2 2 
x_=m + A) 2 O., Xx, =m, + 4 | A o 
2 2 i 2 2 


(3.15) 
h 1 


1 E 
2 2 2 
(2) 

= 2 pz | 


W przypadku kiedy współczynnik skośności y =0 (rozkład zmiennej X jest 


symetryczny) rozwiązanie upraszcza się i ma następującą postać: 


1 
u Mig EE: (3.16) 


Jeżeli do równania (3.11) wstawimy uzyskane wielkości, oraz przyjmiemy N=2 to 
otrzymamy : 
E|Y* |= g" (x_)P +g'(x,)P,, (3.17) 
a zatem dwa pierwsze momenty probabilistyczne mają następującą postać: 
m, =E|Y'|=y_P +y,P, 
2 2 2 3 2 2 DAR 
o, = E|Y |-m =y P +y ml, 
gdzie y, =g(x,) oraz y =g(x_). 
W przypadku gdy Y jest funkcją n nieskorelowanych, ciągłych zmiennych: 
X,,X,,..,X, (każda o zerowym współczynniku skośności), możliwe jest zastąpienie 
każdej z nich zmienną dyskretną. Prawdopodobieństwo (wagę) wyznacza się dla dwóch 


wartości: x_ 1 x, (mniejszej i większej od wartości oczekiwanej o wielkość odchylenia 


standardowego) dla każdej ze zmiennych - łącznie do wyznaczenia potrzebne jest 2” 
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punktów. Każdemu z tych punktów przypisane jest prawdopodobieństwo, które oblicza 


się następująco [Christian i inni, 1999]: 


i TERN 2" 2 > (s,)(5,)rz, (3.19) 


-1 dla x, =m, -o, 
gdzie : s,= | | 


„a ryy Oznacza współczynnik wzajemnej 
+1 dla x, =m, +o, l 


korelacji pomiędzy zmiennymi X, a X. 
W przedstawionym podejściu momenty probabilistyczne losowych funkcji w 
przybliżeniu stanowią sumę wszystkich możliwych kombinacji iloczynów wartości 
funkcji y, i przypisanej jej wagi (prawdopodobieństwa): 

E|Y* ]=>'(y,)P (3.20) 
W przypadku dwóch (n=2) zmiennych: X,,X, wartości prawdopodobieństw 


wyznaczone zgodnie z formułą 3.19 wynoszą: 


F sss) = [1+ (s) Cs) rn ] (3.21) 


Prawdopodobieństwa wyznacza się dla lokalizacji punktów pomniejszonych (—) i 
powiększonych (+) o wartość odchylenia standardowego względem wartości średniej, 
dlatego też w przypadku dwóch zmiennych konieczne jest wyznaczenie 4 wartości 
prawdopodobieństwa : P_ dla punktu o współrzędnych (m, —0, ,m, -0,„), P_ dla 
(m, +0,,m, -0,) , P, dla (m, -0,,m, +0, ) oraz P., dla (m, +0, „m, +0,). 


Zgodnie z formułą 3.21 prawdopodobieństwa wynoszą: 
1 
P.=P_=—|l+r.|, P_=P.=—l1-r. 3.22 
a er en [eral (3.22) 


W przypadku, jeśli zmienne X, i X, nie są skorelowane (ry = 0), 
prawdopodobienstwa dla wszystkich punktów sq takie same i równe 0,25. 
Wartość oczekiwaną oraz wariancję zmiennej Y, wyznaczone na podstawie formuły 
3.20, można wyrazić następująco: 

m, = Yy Pa +y TYT RY 


(3.23) 
Žo ad 2 2 2 2 
6, ia b + yy + W + y_P_ m m, 


gdzie Vax = E(X X) > Y- = BOM) Ya = g(%_,X,,) EJ Y = 9(X,_,X,_). 


OSLWIELUZY.D 
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Lokalizacja punktów oraz przypisane im wagi dla przypadków: dwóch i trzech 


zmiennych są przedstawione na rysunku 3.1. 


P.. 
A X% AX, A + 
(1-1x,.x2)/4 (1+rx,.x)/4 
? — H Y P+++ 
P, + 
P... 
My, È a 
N 
Riak xA e-~« P, s . P+. 
| 20%, J (H-Fxqx2)/4 20y, Y 
[s 1 
T > b d > x 
Mx, 
a) zmienne: X,,X2 X, b) zmienne: X;,,X,,X3 


Rysunek 3.1 Lokalizacja punktów oraz przypisane im wagi dla przypadków: dwóch (a) i trzech (b) 
zmiennych losowych 
W przedstawionym standardowym podejściu Rosenblueth'a liczba punktów 
potrzebnych do oszacowania wynosi 2”. W przypadku dużej liczby n zmiennych 
procedura PEM przestaje stanowić efektywną alternatywę dla innych metod ze względu 


na pracochłonność obliczeń. 


Rosenblueth zaproponował, aby w przypadku gdy zmienne nie są skorelowane 
obliczenia wykonać dla (2n+1) punktów (wyznacza się trzy wartości funkcji Y dla 
każdej ze zmiennych) [Nowak i inni, 2000]: 

Yo =y(my,,,my,,...,my, ) 
y; =g(m, My, ,...m, +0x,..,m, | (3.24) 
Y; Se | Mys -0x,--M, ) 


Na ich podstawie szacuje się wartość średnią Y oraz wariancję V, wykorzystując 


następujące procedury: 
A = RE me (3.25) 
Yo YoYo Yo 
1+7; =(1+V; |(1+V; )...(1+V7 ) (3.26) 
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W powyższych równaniach Y, oznacza wartość średnią wyznaczoną dla wszystkich 
zmiennych poza i-tą, która przyjęta jest jako stała i równa jej wartości oczekiwanej, 
natomiast V, jest współczynnikiem wariancji wyznaczonym w ten sposób, że i-ta 
zmienna przyjęta jest jako losowa a wszystkie poza nią są stałe i równe ich wartościom 
oczekiwanym. 

Redukcję liczby szacowanych punktów z 2” do 2n, umożliwia zastosowanie 
algorytmu zaproponowanego w pracy [Hong, 1998]. Metoda ta może być zastosowana 
tylko w przypadku, gdy zmienne nie są skorelowane, natomiast ich rozkłady może 


charakteryzować brak symetrii. Dla każdej ze zmiennych wyznacza się wielkości: 


(3.27) 
y; = £ (my, > My, SEE, veo My, ) 
na podstawie których szacuje się punkty oraz przypisane im prawdopodobienstwa: 
y | 
x, z m, +0, X; + n | Xi | (3.28) 
X; “| 2 2 
p= Ilya 1 (3.29) 


Zaletą metody estymacji punktowej (PEM) jest niewielka, w porównaniu np. do 
metody Monte Carlo, liczba realizacji potrzebnych do otrzymania estymatorów 
poszukiwanych wielkości. To, jak duża liczba realizacji jest wymagana zależy od liczby 
zmiennych losowych. W przypadku występowania większej liczby zmiennych możliwe 
jest zastosowanie modyfikacji pierwotnego podejścia Rosenblueth’a umożliwiających 
zmniejszenie liczby próbek z 2” do 2n lub 2n+1. Ponadto metoda PEM nie wymaga 
pełnej wiedzy o rozkładach zmiennych- wystarczy znajomość kilku pierwszych 
momentów statystycznych (w podejściu Honga). 

Wadą metody jest brak możliwości weryfikacji wyników obliczeń. Wykonywanie 
analizy jakiegokolwiek problemu inżynierskiego jedynie za pomocą metody PEM wiąże 
się z dużym prawdopodobieństwem uzyskania błędnych rozwiązań. Konieczna jest 


dodatkowa weryfikacja dokonana innymi metodami. 
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3.4 Techniki redukcji wariancji 
3.4.1 Prébkowanie warstwowe (Stratified Sampling SS) 


W przypadku techniki próbkowania warstwowego przestrzeń zmiennych 
losowych Q podzielona jest na pewną ilość rozłącznych podprzestrzeni Q, (/=1,2,..., L) 
spełniających warunek: 

0-Ua ,QNO,=O, (i+ j; i= lene l=1,2,...,L) (3.30) 
I=1 
W powyższej formule L oznacza liczbę wydzielonych podprzestrzeni. 
Rozmiar i kształt warstw przestrzeni Q. mogą być dobrane na 3 sposoby [Shields i inni, 
2015]: 

1.  Symetrically Balanced Stratified Design (SBSD): przypadek, gdy występuje 
równe prawdopodobieństwo znalezienia się w danej warstwie oraz przestrzeń 
zdarzeń podzielona została na warstwy o równych wymiarach (warstwy są n- 
wymiarowymi hiperszescianami). 

2. | Asymmetrically Balanced Stratified Design (ABSD): przypadek, gdy występuje 
równe prawdopodobieństwo znalezienia się w danej warstwie natomiast 
warstwy maja różne wymiary (warstwy sq  n-wymiarowymi 
hiperprostopadłościanami). 

3. Unbalanced Stratified Design (UBSD): przypadek, w którym występują różne 
prawdopodobieństwa znalezienia się w danej warstwie oraz wymiary warstw nie 
są sobie równe. 

Przykładowy podział przestrzeni zdarzeń 62 na warstwy dla przypadku zagadnienia 


dwuwymiarowego) przedstawiono na rysunku 3.2. 


Rysunek 3.2 Przykładowe sposoby wyodrębnienia warstw z dwuwymiarowej przestrzeni Q 
zmiennych [Zhang i inni, 2010] 
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Jeżeli przez p, (I=1,2,...,L) oznaczone zostało prawdopodobieństwo znalezienia się w 
warstwie Q; , natomiast p, wyraża prawdopodobieństwo awarii realizacji pochodzącej 


z tej warstwy, to całkowite prawdopodobieństwo awarii dla /-tej podprzestrzeni można 
wyrazić jako: 
pi= | FO, (1=123,.,1) (3.31) 
Q; 


Dla warstwy Q, funkcja gęstości rozkładu prawdopodobieństwa ma postać: 


0 dla x€Q, 
= FH 12 Beale 3.32 
Sf, (x) f(x) dla seam i ee aa ) ( ) 
Pi 


W przypadku tak zdefiniowanej funkcji /, (x), prawdopodobieństwo awarii realizacji 


pochodzącej z podprzestrzeni Q, szacuje się jako: 


— Ni lo, (x,,)-F(x,,) s 
By > fils," (1 =1,2,3,...,L) (3.33) 


W powyzszej formule x,, oznacza j-tą wygenerowaną (zgodnie z funkcją /,(x)) 


próbkę pochodzącą z wektora realizacji, z podprzestrzeni Q,, natomiast N, oznacza 


całkowitą liczbę realizacji w obrębie /-tej warstwy. 


Całkowite prawdopodobieństwo awarii szacowane dla podprzestrzeni wyraża formuła: 


L LIN I x, |- f {x 
BP = Diy > SPORE , (1=1,2,3,...,Ł) (3.34) 
Błąd powyższego estymatora prawdopodobieństwa awarii wyraża się poprzez jego 
wariancję, analogicznie jak miało to miejsce w przypadku bezpośredniej metody Monte 
Carlo (formuła 3.6). 

Należy podkreślić, że losowanie warstwowe jest techniką umożliwiającą w 
znacznym stopniu zredukowanie czasu obliczeń numerycznych, jednak aby tak się stało 


proces podziału przestrzeni na warstwy musi być właściwie przeprowadzony. 


Pobrano z mostwiedzy.pl 
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3.4.2 Technika próbkowania hipersześcianem łacińskim (Latin Hypercube 
Sampling (LHS)) 


W metodzie próbkowania hipersześcianem łacińskim, podobnie jak w 
próbkowaniu warstwowym, przestrzeń zmiennych dzielona jest na rozłączne 


podprzestrzenie Q, (/=1,2,..., L), jednak w przypadku próbkowania LHS generowanie 


próbek z utworzonych podprzestrzeni podlega pewnemu ograniczeniu. Jeżeli przestrzeń 
jest n- wymiarowa oraz dla każdego i-tego wymiaru następuje podział na / 
podprzedziałów, to z każdej utworzonej w ten sposób podprzestrzeni wylosowana może 
być pojedyncza próbka pozyskana na zasadzie kombinatoryki pomiędzy wszystkimi n 
wymiarami wektora zmiennych [McKay i inni, 1979]. 

Określenie „hipersześcian łaciński” jest rozszerzeniem pojęcia kwadratu 
łacińskiego Eulera, czyli /-poziomowej dwuwymiarowej tablicy wypełnionej liczbą Z 
wartości usytuowanych w taki sposób, że żaden wiersz ani kolumna nie zawiera dwóch 


tych samych wartości. Przykłady kwadratów łacińskich przedstawiono na rys. 3.3. 


Rysunek 3.3 Przykłady kwadratów łacińskich [źródło: www.bing.com/images, 2018] 


W metodzie LHS, w przypadku n- wymiarowej przestrzeni realizacji, w której każdy z 
wymiarów podzielony został na / równie prawdopodobnych podprzedziałów, liczbę 


możliwych kombinacji próbkowania wyznaczyć można z następującej formuły: 


L -|fit- i) = (1!)"" (3.35) 


j=0 


= 
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Próbki z poszczególnych podprzedziałów każdego z n wymiarów przestrzeni realizacji 
można dobrać jako odwrotność dystrybuanty rozkładu prawdopodobieństwa kolejnych 
składowych wektora X, (rys.3.4). Związek ten opisuje następujący wzór: 


_1( k-0,5 
al ? | (3.36) 


W powyższej formule x, oznacza k- tą próbkę i-tej składowej wektora X, , natomiast L 


jest całkowitą liczbą próbek możliwych do uzyskania dla i-tego wymiaru wektora 


realizacji zmiennych [Helton i inni, 2003]. 


F 


Xi 


>x 


ik 


Rysunek 3.4 Schemat wyboru k- tej próbki i- tej zmiennej wektora X, e ) na podstawie 


dystrybuanty jej rozkładu prawdopodobieństwa 


W przypadku liczby Z, próbek wektora zmiennych losowych X , wyznaczonych na 


drodze próbkowania metodą LHS, estymator prawdopodobieństwa awarii szacuje się 


następująco: 


(3.37) 


gdzie n( g(x, ) < 0) oznacza liczbę symulacji, dla których konstrukcja uległa awarii 


natomiast M, oznacza całkowitą liczbę symulacji (N, < L, ). 


Tak jak w przypadku innych metod probabilistycznych radykalnie 
ograniczających liczbę próbek obliczenia wykonane za pomocą LHS nie dają możliwości 
weryfikacji wyników i powinny być sprawdzane innymi dostępnymi algorytmami 


szacującymi niezawodność oraz inne parametry opisujące konstrukcję. 
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3.4.3 Technika próbkowania ukierunkowanego (Targeted Random Sampling 
TRS) 


Jednym z wyzwań we współczesnej analizie inżynierskiej jest ocena 
prawdopodobieństwa awarii systemów, których funkcja stanu granicznego jest złożona 
(nieliniowa, nieciągła, niemonotoniczna). Przez ostatnie dekady opracowano szereg 
metod umożliwiających oszacowanie niezawodności takich układów. Ważnym 
czynnikiem przeważającym o ich stosowalności jest łatwe zaaplikowanie danej metody 
w programie komputerowym. 

W poniższym rozdziale zaprezentowano podejście tzw. ukierunkowanego 
losowego pobierania próbek (TRS), zaproponowane przez [Shields i inni 2015 a]. Metoda 
ta umożliwia oszacowanie prawdopodobieństwa awarii w przypadku problemów ze 
złożonym stanem granicznym. Metoda wykorzystuje tzw. doskonaloną warstwową 
koncepcję pobierania próbek RSS (Refined Stratified Sampling) opisanej w [Shields i 
inni, 2015 b], której szczegółowy algorytm przedstawiono w dalszej części pracy. 

Technika próbkowania ukierunkowanego umożliwia wyodrębnienie spośród 
początkowo przyjętej populacji tych próbek, które skoncentrowane są w bezpośrednim 
sąsiedztwie stanu granicznego. Pozwala to na określenie obszaru awarii Q, a co za tym 
idzie oszacowanie prawdopodobieństwa p, , przy użyciu niewielkiej liczby symulacji. 

Wydzielenie docelowej populacji próbek w technice TRS odbywa się poprzez 


kolejne podziały wstępnie założonych, rozłącznych warstw. Rozpatrując dwuwymiarową 


przestrzeń zmiennych losowych X=[X,,X,], algorytm techniki TRS można 


przedstawić następująco: 
1. Wstępny podział przestrzeni zdarzeń (2 na rozłączne warstwy (rysunek (3.6a)) 


wykorzystujący technikę próbkowania warstwowego oraz klasyfikacja zdarzeń z 
poszczególnych warstw jako awaryjne (f x;) lub bezawaryjne (x, ). Początkowy 
podział musi być przeprowadzony w taki sposób, aby istniała co najmniej jedna 
warstwa zawierająca zdarzenie awaryjne f Q. 


2. Identyfikacja podprzestrzeni próbkowania celowego: poszukiwanie sąsiadujących 


ze sobą (w dowolnym ortogonalnym kierunku) par warstw, z których jedna zawiera 


zdarzenie awaryjne / Q, , druga zaś zdarzenie bezawaryjne °Q j» Jeżeli obszar, jaki 
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zajmuje dana para warstw oznaczymy jako ©; = f oU Q j » to całkowity zasięg 


próbkowania celowego wynosi I = U ©; . Na rysunku 3.6b obszar I’ uwydatniono 
ij 


poprzez zakreskowanie. 


Podziaął warstw zgodny z następującą procedurą: spośród wszystkich 
wyznaczonych par warstw ©, wybrana zostaje taka ©”” =©,,, dla której 
zachodzi: 
P| 6% |-P[e,]=P| "o,lJ'Q, |zP[e,|=P|"QlJ'e,| 639 
a następnie konstruowany jest wektor łączący punkty próbek danej pary warstw 
By =" Xx, "X; (3.39) 
oraz wyznaczone są jego rzuty na dwa ortogonalne kierunki: 


=é,-*x,, (3.40) 


l 


Rzut wektora * x,, 0 większej długości wyznacza kierunek podziału warstwy. 
Odbywa się on poprzez poprowadzenie prostej prostopadłej do dłuższego z rzutów, 


przechodzącej przez punkt przecięcia odcinka * x,, z funkcją stanu granicznego 


G(x)=0 (rysunek 3.5): 


Kewl, =0( x) (3.41) 
Warstwa zawierająca wyznaczony punkt przecięcia zostaje podzielona a obszar 
próbkowania zwiększa się o kolejną podprzestrzeń a zatem również o kolejny 
wynik symulacji. Proces ten jest powtarzany tak długo, aż wyznaczona zostanie 
funkcja stanu granicznego i obszar awarii. Koncepcję podziału warstw 


przedstawiono na rysunku 3.6 g-i. 
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Rysunek 3.5 Ilustracja długości odcinka x,,,, wyznaczającej lokalizację podziału warstwy (na 


podstawie [Shields i Sundar, 2015] 


4. Oszacowanie prawdopodobieństwa awarii : 
Nf 
pr => P £9, | (3.42) 
i=l 


W powyższej formule N, oznacza całkowitą liczbę zdarzeń awaryjnych 


pochodzącą z rozpatrywanego zbioru próbek. 


Podstawą właściwie przeprowadzonej techniki TRS jest wstępny podział 

przestrzeni próbkowania na warstwy oraz poprawne sklasyfikowanie pochodzących z 
nich zdarzeń. 
Awaria konstrukcji nie zawsze jest zjawiskiem łatwo przewidywalnym- ważne jest, ażeby 
wyjściowe uwarstwienie zostało przeprowadzone w taki sposób, aby uwzględniało 
możliwość wystąpienia zdarzeń będących kombinacją różnych wartości zmiennych 
(niekoniecznie wartości ekstremalnych, pochodzących z ogonów ich rozkładów 
prawdopodobieństwa) [Shields, 2014]. 

Zastosowanie techniki próbkowania ukierunkowanego, nawet dla przypadków 
nieliniowej funkcji stanu granicznego prowadzi do szybkiego wzrostu dokładności 
szacowanego prawdopodobieństwa awarii przy stosunkowo niewielkiej liczbie próbek. 
Istotne ograniczenie stanowi natomiast wymiar wektora zmiennych losowych- technika 


ta sprawdza się w problemach o małej liczbie zmiennych. 
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Rysunek 3.6 Technika próbkowania ukierunkowanego dla przypadku dwuwymiarowego: ogólna 

koncepcja metody(a-f), szczegółowa prezentacja podziału warstw (g-i) [Shields i Sundar, 2015] 
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3.4.4 Technika próbkowania (Refined Stratified Sampling RSS) 


Zastosowanie metody RSS (Refined Stratified Sampling) ma na celu zagęszczenia 
próbek w przestrzeni zmiennych Q. 
Rozpatrywany jest wektor zmiennych losowych x zawierający n nieskorelowanych 
zmiennych projektowych zdefiniowanych w przestrzeni (2. Biorąc pod uwagę zbiór N 
próbek podzielonych na M warstw (M=N) metodologia RSS przebiega w następujący 
sposób [Shields i inni, 2015 b]: 


1. Wybór warstwy, która zostanie podzielona -Q ; spełnione muszą być następujące 
warunki: 
= jeżeli istnieje jedna warstwa Q, taka, że w,>w, V jek to może ona zostać 
podzielona, 


= jeżeli istnieje N, warstw Q,, k =1,2,...,N,, których waga w, = max (w A to jedna 
J 


z nich zostaje wybrana do podziału w sposób losowy z równym 
prawdopodobieństwem równym P[D(Q,)]=1/N, 
2. Podział wybranej warstwy ©, na pół zgodnie z następującymi wytycznymi: 

= wyznaczenie długości krawędzi warstwy przed podzieleniem we wszystkich jej 
kierunkach: A, = é% -ć ,i=1,2,...,n, 

" wybór najdłuższego boku warstwy przed podziałem 4, = max (Ax); 

= podział warstwy na pół w kierunku prostopadłym do najdłuższego wymiaru 
warstwy -4, ; jeżeli istnieje M, krawędzi warstwy, których długość równa jest 


A, =A,, i=1,2,...,N., to wybór jednej z nich następuje w sposób losowy z 


prawdopodobieństwem równym P[D(i “\}=1/ N. 
3. Ponowne wyznaczenie wag warstw po podziale. 
4. Powtórzenie procedury opisanej w punktach 1-3 aż do uzyskania pożądanej liczby 


próbek 


Przykład podziału metodą RSS dla przypadku dwóch zmiennych i podziału 
rozpoczynającego się od jednej warstwy (M=N=l1) przedstawiono na rys. 3.7. Schemat 


przedstawia przebieg czterech pierwszych podziałów. 
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6. podział warstwy o większych 
wymiarach Q,.. 
podział przestrzeni zmiennych 
na warstwy Q,;k=1, 2, 3, 4 


2. losowy podział 
przestrzeni zmiennych 
na warstwy Q,.;k=1,2 


x lw 


1. wylosowanie próbki 
z przestrzeni zmiennych 


3. losowy podział 
przestrzeni zmiennych 
na warstwy Q,.;k=1, 2, 3 


8. losowy podział 
przestrzeni zmiennych 
na warstwy Q,.; 
k=1, 2, 3,4,5 


====uz 
`~ 
`~ 
`~ 


A 


3. wylosowanie próbki 
w pustej warstwie 
oraz przypisanie nowych wag 
w,=P(Q,) 


7. wylosowanie próbki 
w pustej warstwie 
oraz przypisanie nowych wag 
w,=P(Q,) 


5. wylosowanie próbki bej" we 

w pustej warstwie I 
oraz przypisanie nowych wag |. | 
w=P(Q,) Ll 


9. wylosowanie próbki 
w pustej warstwie 
oraz przypisanie nowych wag 
w,=P(Q,) 


Rysunek 3.7 Metoda RSS (Refined Stratified Sampling)- schemat kolejnego podziału przestrzeni na warstwy 


|d'‘Azpaimjsow z ouejqod A ZOJIM LI SOW RZY 
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Rozdział 4. Analiza wrażliwości konstrukcji 


4.1 Wprowadzenie 


Analizą wrażliwości nazywamy badanie wpływu zmiany parametrów 
wejściowych modelu matematycznego na zmienność jego wyników. Pojęcie analizy 
wrażliwości często mylone jest z pojęciem analizy niepewności modelu. Zadaniem 
pierwszej jest badanie w jaki sposób niepewność odpowiedzi układu można przypisać 
różnym jej źródłom zawartym w danych wejściowych, natomiast druga zajmuje się 
charakterystyką niepewności odpowiedzi. Bardzo często analiza niepewności stanowi 
wstęp do analizy wrażliwości (najpierw wyznacza się charakterystyki zmiennej 
wyjściowej, po czym próbuje się podzielić źródła jej niepewności według wagi). 

Analiza wrażliwości konstrukcji możliwa jest do przeprowadzenia zarówno 
metodami deterministycznymi jak również probabilistycznymi. Obecnie powszechnie 
stosowane w projektowaniu są obliczenia deterministyczne. Przeprowadza się je jako 
analizy parametryczne, wykorzystuje się również metody wariancji lub określa efekt 
wpływu pojedynczych parametrów na zmienne wyjściowe. Tego typu podejście nie daje 
jednak pełnej informacji o wpływie imperfekcji geometrycznych 1 materiałowych 
elementów, a także niepewności opisu obciążenia na mechaniczną odpowiedź 
konstrukcji. Poszukiwany jest szerszy — globalny opis wrażliwości uwzględniający 
możliwą zmienność wszystkich parametrów definiujących konstrukcję [Chan i inni, 
1997]. Wydaje się więc istotne włączenie metod probabilistycznych a także 
stochastycznych do procesu bezpiecznego i ekonomicznego projektowania. Tego typu 
analizę można określić jako wrażliwość probabilistyczną lub niezawodnościową, 
poszukującą wpływu zmienności losowej każdej ze zmiennych podstawowych zadania 
na losową niezawodność elementu/ukladu. Takie sformułowanie umożliwia 
przeprowadzenie zarówno analizy poziomu II, zdefiniowanej na poziomie wartości 
oczekiwanych i odchyleń standardowych, jak również poziomu III, opisanej uzyskanymi 
na drodze numerycznej histogramami. 

W metodach poziomu II funkcja stanu granicznego (liniowa lub zlinearyzowana) 
zawiera wpływy liniowe poszczególnych zmiennych podstawowych. Współczynniki 
wrażliwości (sensitivity factors) utożsamiane są zazwyczaj ze współczynnikami 
związanymi z liniowymi zmiennymi podstawowymi. Można je również zdefiniować jako 


pochodne cząstkowe funkcji stanu względem danej zmiennej w punkcie wartości 
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średnich lub w punkcie projektowym [Nowak, 2000]. Ocena wrażliwości 
przeprowadzana na podstawie uzyskanych histogramów funkcji stanu granicznego 
(poziom III), uzależnionych od poszczególnych zmiennych projektowych, nie jest łatwa. 
Wynikowy histogram zazwyczaj nie ma charakteru funkcji Gaussa, a każda próba 
aproksymacji nie jest jednoznaczna i zmienia rozwiązanie problemu. Ten kierunek 


badania wrażliwości wymaga dalszych analiz. 


Analiza wrażliwości wykonywana w procesie tworzenia modelu matematycznego 


umożliwia m.in. [Hamby, 1994], [Saltelli i inni, 2000]: 


u wyeliminowanie spośród często licznej grupy parametrów projektowych tych, 
które nie wpływają znacząco na wynik. Pominięcie tego kroku prowadzi do 
zwiększenia czasu obliczeniowego i nieuzasadnionego wydłużenia obliczeń; 

= wyodrębnienie spośród wszystkich parametrów tych, których zmiana w 
największym stopniu wpływa na wynik; 

=" określenia wpływu korelacji pomiędzy poszczególnymi parametrami; 

« sprawdzenie czy stworzony model matematyczny jest wiarygodny np. poprzez 
skontrolowanie czy na skutek zmiany przypuszczalnie nieistotnych parametrów nie 
występują silne wahnięcia wyniku lub zweryfikowanie czy zakres tej zmiany jest 


racjonalny. 


Celem analizy wrażliwości jest zatem rozpoznanie mechanicznego zachowania się 


konstrukcji i możliwe uproszczenie modelu matematycznego. 
4.2 Klasyfikacja metod analizy wrażliwości 


Obecnie istnieje obszerne spektrum metod służących analizie wrażliwości, 
wykorzystujących rozwiązania analityczne bądź różnego rodzaju techniki 
probabilistyczne. W literaturze spotkać się można z różnymi klasyfikacjami metod 
analizy wrażliwości, w zależności od charakterystyki rozpatrywanego problemu. W pracy 


[Hamby, 1994] wyróżniono trzy grupy metod: 


1. te, które oceniają wynik w przypadku zmiany jednego parametru na raz, między 


innymi : analiza różniczkowa DA (ang. Differential Analysis), technika OAT (One- 
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At-a-Time), metody: FA (Factorial Analysis), IF (Importance Factors) oraz SI 
(Sensitivity Index); 

te, które opierają się na generacji próbek wektorów zmiennych wejściowych i 
identyfikacji powiązanych z nimi odpowiedzi układu, m. in.: metody analizy 
regresji, stosowanie wykresów rozproszenia, wyznaczanie współczynników: 
korelacji Pearsona lub korelacji cząstkowej; 

te, które wymagają odpowiedniego podziału wektora zmiennych wejściowych na 
podstawie uzyskanego wektora odpowiedzi, m in. test Kolomogrowa-Smirnowa, 


test Cramera-von Misesa, test Manna-Whitneya. 


Inną klasyfikację metod znaleźć można w [Frey i inni, 2003]. Rozróżnia ona: 


Matematyczne metody analizy wrażliwości, umożliwiające analizę modeli 
deterministycznych, polegające na przeprowadzaniu serii badań numerycznych 
uwzględniających zmienność pojedynczego czynnika (w jego możliwym zakresie) 
i badaniu wpływu tej zmiany na wynikową odpowiedź układu. W metodach 
matematycznych wrażliwość jest najczęściej wyrażona poprzez wartość względną 
zmiany odpowiedzi układu. Głównym celem stosowania tej grupy metod jest 
szacowanie wpływu zakresu zmienności parametrów wejściowych na odpowiedź 
układu oraz identyfikacja parametrów, których zmiana ma największe 
oddziaływanie na końcowy rezultat [Morgan i inni, 1990]. Do metod 
matematycznych zaliczyć można między innymi: NRSA (Nominal Range 
Sensitivity Analysis) oraz DSA (Differential Sensitivity Analysis). 
Probabilistyczne (lub statystyczne) metody oceny wrażliwości wykorzystujące 
zbiór wyników symulacji, przeprowadzonych na podstawie zmiennych 
wejściowych, charakteryzujących się założonymi rozkładami 
prawdopodobieństwa. W zależności od przyjętej metody badana jest odpowiedź 
układu na zmianę jednej bądź kilku zmiennych wejścia na raz. Metody 
probabilistyczne umożliwiają ilościową ocenę wpływu jednoczesnej zmiany kilku 
parametrów na odpowiedź układu. Wśród metod probabilistycznych wymienić 
można m.in.: analizę regresji liniowej LRA (Linear Regression Analysis), analizę 
wariancji ANOVA (ANalysis Of VAriance) metodę powierzchni odpowiedzi RSM 
(Responce Surface Method), metodę FAST (Fourier Amplitude Sensitivity Test). 
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Graficzne metody oceny wrażliwości wykorzystywane często jako uzupełnienie 
metod matematycznych i statystycznych, umożliwiające łatwiejszą interpretację 
uzyskanych wyników. Użycie graficznej reprezentacji danych (w postaci np.: 
wykresów, map lub powierzchni) może stanowić metodę przesiewową przed dalszą 
analizą modelu jak również służyć obrazowaniu zależności pomiędzy zmianą 
poszczególnych parametrów a odpowiedzią układu [Frey i inni, 2003]. Przykładem 
wykorzystania metody graficznej jest stosowanie histogramów, wykresów 


rozproszenia (scatterplots) lub wykresów pudełkowych (boxplot). 


Trzecim rodzajem klasyfikacji metod analizy wrażliwości przedstawionym między 


innymi w [Heiselberg i inni, 2009 ] oraz [Saltelli i inni, 2004] jest podział na: 


1. 


Lokalne metody analizy wrażliwości, wykorzystujące często na podejście OAT 
(One-At-A-Time), skupiające się na badaniu wpływu odchyleń pojedynczego 
parametru od jego wartości nominalnej (może to być np. wartość Średnia) na 
odpowiedź układu przy zachowaniu stałych wartości pozostałych czynników. 
Metody lokalne, w przypadku gdy model jest nieliniowy lub jeśli różnice 
zaobserwowane w odpowiedzi układu powstałe na skutek modyfikacji 
poszczególnych parametrów znacząco różnią się od siebie, powinny być zastąpione 
metodami globalnymi [Cukier i inni, 1973]. Przykładem metody lokalnej jest 
podejście DA (Differential Analysis) wykorzystujące obliczenie (lub oszacowanie) 
pochodnych cząstkowych funkcji stanu modelu w konkretnym punkcie przestrzeni 
zmiennych [Bertrand i inni, 2015]. 

Globalne metody analizy wrażliwości, których celem jest badanie, w jaki sposób 
niepewność występująca w danych wyjściowych modelu może być pogrupowana 
według źródeł niepewności uwzględnionych w danych wejściowych [Saltelli, 
2002]. Metody te umożliwiają analizę zachowania się odpowiedzi modelu na 
skutek modyfikacji dowolnej liczby jego parametrów wejściowych w całej 
przestrzeni ich zmienności. W podejściu opartym na globalnej analizie wrażliwości 
wpływ danego parametru początkowego na odpowiedź modelu jest oceniany przy 
jednoczesnej zmianie wszystkich innych czynników wejściowych [Nguyen i inni, 
2015]. Do grupy metod globalnej analizy wrażliwości zaliczyć można między 
innymi: analizę FAST (Fourier Amplitude Sensitivity Test), techniki mające za 


podstawę m. in. metodę Monte Carlo lub metodę powierzchni odpowiedzi. 
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Niezależnie od przyjętej metody, proces globalnej analizy wrażliwości można sprowadzić 

do następujących etapów: 

Etap 1: definicja modelu oraz ustalenie, które spośród zmiennych wejściowych zostaną 
poddane analizie; 

Etap 2: dobór parametrów probabilistycznych lub zakresów zmienności każdej ze 
zmiennych wejściowych. Podstawą do ustaleń może być m.in.: wiedza 
zaczerpnięta z literatury, dane doświadczalne lub opinia ekspertów [Saltelli i 
inni, 2008]; 

Etap 3: generacja wektora (lub macierzy) zmiennych wejściowych; 

Etap 4 obliczenia numeryczne modelu; 

Etap 5: ocena wpływu zmiany rozpatrywanego parametru na podstawie danych 
wyjściowych modelu, wyciągnięcie wniosków i decyzja o ewentualnym 


rozpoczęciu nowej iteracji analizy. 
4.3 Wybrane metody analizy wrażliwości 
4.3.1 Projektowanie przesiewowe (Screening Design) 


Terminem ,,projektowanie przesiewowe” określa się eksperymenty numeryczne, 
wykonywane na wstępie właściwej analizy modelu, mające na celu wyodrębnienie 
(spośród często licznego zbioru) tych czynników, które w największym stopniu wpływają 
na odpowiedź układu. Szeroki opis metod badań przesiewowych można znaleźć m.in. w 
pracy [Woods i inni, 2015]. 

Algorytmy obliczeniowe metod przesiewowych powinny być zorganizowane w 
taki sposób aby umożliwiały analizę dużej liczby zmiennych, a więc powinny być 
obliczeniowo i czasowo wydajne. Często metody te dostarczają jedynie informacji o tym, 
w jaki sposób parametry projektowe szeregują się pod względem ważności, natomiast nie 
udzielają odpowiedzi na pytanie o ile dany czynnik jest bardziej istotny niż pozostałe 
[Saltelli i inni, 2000]. 

W badaniach przesiewowych przeważnie wykorzystują podejście OAT. W tego 
typu analizie eksperymenty wykonuje się dla ustalonych trzech poziomów rozpatrywanej 
zmiennej: wartości ekstremalnych oraz wartości środkowej (znajdującej się „w połowie 
drogi” pomiędzy wartościami ekstremalnymi). Eksperyment przeprowadzony z 


wykorzystaniem wartości środkowych zmiennych zwany jest eksperymentem 
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kontrolnym. Ocena parametrów, na zmianę których model jest wrażliwy polega na 
wyznaczeniu i porównaniu ze sobą różnic między wynikami uzyskanymi dla 
ekstremalnych poziomów zmiennych z tymi uzyskanymi w eksperymencie kontrolnym 
[Saltelli i inni,2000]. Istotnym ograniczeniem podejścia OAT jest fakt, że w tego typu 
analizie szacuje się jedynie efekt główny danej zmiennej - wpływ potencjalnie 
występującej między dwoma lub więcej czynnikami interakcji jest całkowicie pominięty. 

Inną metodą badań przesiewowych jest zastosowanie tzw. eksperymentowania 
czynnikowego (Factorial Experimentation), które umożliwia wyznaczenie odpowiedzi 
układu dla wszystkich możliwych kombinacji przyjętych poziomów zmiennych. 
Wyróżnia się badanie czynnikowe pełne i frakcyjne. W przypadku pełnego eksperymentu 
czynnikowego dla przyjętych k poziomów każdej z N zmiennych losowych całkowita 
liczba eksperymentów wynosi N*. W przypadku eksperymentowania frakcyjnego 
wykorzystuje się tylko wybrane poziomy poszczególnych zmiennych. Dobór punktów 
ustalonych w badaniu czynnikowym: pełnym i frakcyjnym dla przypadku 


trójwymiarowego zobrazowano na rys. 4.1. 


AX. . o . X . 


va 5 > 


X, X, 


(a) (b) 


Rysunek 4.1 Przykładowy wybór punktów eksperymentowania czynnikowego pełnego (a) i 
frakcyjnego (b) — przypadek trójwymiarowy [źródło: 
https://www.slideshare.net/maheshmarathe4/factorial-design-69296752, 2018] 


Kolejną gałęzią metod przydatnych w badaniach przesiewowych są metody 
graficzne, obrazujące wrażliwość modelu na zmianę rozpatrywanych parametrów w 
postaci m.in. wykresów, map punktów czy też histogramów. 

Metody graficzne stanowią formę jakościowego określania wrażliwości modelu 


na zmianę rozpatrywanych parametrów wejściowych. Mogą być one stosowane zarówno 


Pobrano z mostwiedzy.pl 
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jako metoda przesiewowa stosowana przed dalszą analizą modelu w celu lepszego 
rozpoznania zależności pomiędzy zmiennymi wejścia i wyjścia, jak również jako 
uzupełnienie metod matematycznych i statystycznych analizy wrażliwości [Frey i inni, 
2003]. Często stosowanym w ramach metod graficznych podejściem jest wykorzystanie 
tzw. chmur rozproszenia punktów. Metoda ta polega na naniesieniu na wykres punktów, 
których pierwszą współrzędną jest wartość wejściowa analizowanej zmiennej, natomiast 
drugą współrzędną jest wynik symulacji (np. bezpośredniej metody Monte Carlo). 
Wykres rozproszenia dostarcza informacji o zakresie wartości dla zmiennych wejścia i 
wyjścia oraz umożliwia identyfikację relacji pomiędzy nimi. Ograniczeniem w 
stosowaniu chmur rozproszenia jest po pierwsze czas ich generacji w przypadku 
skomplikowanych modeli zawierających wiele zmiennych, a po drugie fakt, że ich 
interpretacja oparta jest jedynie na subiektywnym osądzie i może prowadzić do błędów. 
Przykładowy obraz chmur rozproszenia wyników, zaczerpnięty z [Saltelli i inni, 2008] 
przedstawiono na rysunku (4.2a) oraz (4.2b). Zaprezentowane wykresy dotyczą dwóch 
spośród czterech znormalizowanych zmiennych wejściowych Z;, od których zależna jest 


zmienna wyjścia Y. 


40 T T T T T T T 15 a TT U U YU TYT] 
30 - 4 
> 10 
20 - g 5 4 
10 4 80. s ZAM 
WD u u 
>0- 2 -5 4 u 
104 2.104 
-20- >_15 4 
-30 + -20 - 
-40 -25 I ! I I L 
Z, (a) przedziały zmiennej Z, (c) 
40 T_T 15 T Rsi 
30 5 104 
20 - > 5. ET 
„U 
10 7 5 0 4 m 
[> u 
>04 45 -5 - u . 
o 
-10- 8-101 " 
-20 4 $-154 
-30 + -20 - 
-40 -25 I I er ers | | 
Z, (b) przedziały zmiennej Z, (d) 


Rysunek 4.2 Przykładowy obraz chmur rozproszenia zmiennej Y w zależności od wartości 
zmiennych Z; oraz Z? wygenerowanych za pomocą bezpośredniego próbkowania Monte Carlo (rys. 
a i b) oraz wartości średnie zmiennej Y wyznaczone na podstawie zbiorów próbek z wydzielonych 
przedziałów zmiennych Z; i Z2(rys. c i d) [Saltelli i inni, 2008] 
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W przedstawionym przykładzie chmury a i b dotyczą tych samych punktów 
obliczeniowych, różnica w ich kształcie zależy jedynie od tego na jakiej rzutni (Z,Y czy 
Z2Y) zostały one wykreślone. Za najbardziej istotne dla analizowanego problemu uważa 
się te zmienne, których rozproszone punkty charakteryzują się widocznym 
uporządkowaniem oraz największym zakresem zmienności na osi zmiennej wejściowej. 
W tym przypadku zmienną, na którą odpowiedź układu Y jest najbardziej wrażliwa jest 
zmienna Z2. 

Uszeregowanie zmiennych wejściowych pod względem istotności ich zmiany na 
odpowiedź układu odbywać się może zarówno bezpośrednio na podstawie kształtu 
chmury jak również poprzez wyznaczenie linii trendu wartości Średniej zmiennej Y dla 


punktów z wydzielonych podprzedziałów rozpatrywanej zmiennej Z, ( punkty widoczne 


na rys. 4.2 cid). 
4.3.2 Metody badania wrażliwości wykorzystujące analizę wariancji 


Metodę badania wrażliwości wykorzystującą analizę wariancji po raz pierwszy 
zastosował chemik R. I. Cukier we wczesnych latach siedemdziesiątych ubiegłego 
stulecia. Metoda ta, znana pod nazwą FAST (Fourier Amplitude Sensitivity Test) 
umożliwiała określenie głównego efektu rozpatrywanego czynnika na wariancję 
odpowiedzi układu. W latach 90tych rosyjski matematyk I. M. Sobol, rozwijając 
koncepcję, którą zaproponował Cukier, przedstawił procedurę wykorzystującą 
próbkowanie Monte Carlo, umożliwiającą szacowanie miar wrażliwości dla dowolnie 
wybranych grup czynników a przez to uwzględnienie wszelkich interakcji zachodzących 


między czynnikami poddanymi analizie. 
Procedura obliczenia indeksów wrażliwości Sobola metodą Monte Carlo 


Jeżeli model matematyczny opisany jest funkcją f(x), gdzie x=(x%,x,,...,x,) 


S 


jest wektorem zmiennych losowych zdefiniowanych w n -wymiarowej przestrzeni K” , 
takiej że: 

K” ={x|0 <x, <1, i=1,2,..7} (4.1) 
to funkcję f(x) przedstawić można jako następującą sumę: 


FX) 5+205(3)+)) > Í; (x,,%,) ++ fian (X.X, ) (4.2) 


i=l j=i+l 
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Istotną cechą równania (4.2) jest fakt, że wszystkie czynniki w nim zawarte są 


ortogonalne, np. jeżeli (i,,...,i, ) = ( j,,..., j,) to zachodzi: 


E Faun = 0 (4.3) 
a 
Całkowita oraz częściowa wariancja odpowiedzi modelu zdefiniowane są następująco: 
D= I fP (x)dx- fe (4.4) 
1 1 
=] [Ra (0%) dz, dx, (4.5) 
0 0 


gdzie I i, <...<i <n oraz S=l1,2,..,n 
Całkowitą wariancję odpowiedzi modelu przedstawić można w postaci następującej 


sumy: 


n n 


D= p +) > 2, +..+D, n (4.6) 


i=l j=i+1 


Na podstawie powyższego równania wyznacza się indeksy wrażliwości: 
Sissi) = (4.7) 


Suma wszystkich indeksów wrażliwości równa jest 1: 
>,S(i)+V, DY S(i,j)+...+S(1,2,...,n)=1 (4.8) 
i=l i=l j=i+l 


Całki (4.4) oraz 4.5 można wyznaczyć stosując całkowanie Monte Carlo. Oszacowane w 


ten sposób estymatory wielkości f,, Di D, wyrażone są następująco: 


„SEA 4.9 
h yt (ža) (4.9) 
A 1% 2 72 
D- 2) (x, )- fę (4.10) 
m=l 
5 = AOS (x hal F(R )- JE (411) 
m=l 


W równaniach 4.9 — 4.11, N oznacza liczbę generacji, x,, jest punktem obliczeniowym 


m 
r n ed = 
określonym w K”, natomiast x „= M 2 


Indeksy górne ® i ® wskazują, że generowane są dwie macierze o wymiarze N xn 


każda, zawierające wylosowane elementy wektora x). 
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Przedstawiona w pracy [Saltelli i inni, 2008] procedura wyznaczenia indeksów 


wrażliwości Sobola zakłada: 


1. generację dwóch macierzy, każda o wymiarze NXn (n oznacza liczbę zmiennych 
losowych zadania) zawierających wylosowane elementy wektora zmiennych 


(wzory 4.12 i 4.13 ). 


x? x” a” 5, 
(2) (2) (2) (2) 
x M x; n 
A=| .. ie sad ie a 2 (4.12) 
(N-1) (N-1) (N-1) (N-1) 
X Xz i n 
(N) (N) (W) (W) 
X Xz i n J 
a) (1) (1) a) 
Xati Xn? nti A 
(2) (2) (2) (2) 
Xn Xn42 nti Xan 
B=| .. si ssa sę sze "m (4.13) 
W-D)  „(N-I) (N-1) (N-1) 
BZ] Xn+2 R Xati „as Xn 
(N) (N) (N) (N) 
ie Xai Xn42 Xati one Mn 


2. Zdefiniowanie macierzy C,, której wszystkie kolumny poza i-tą pochodzą z 


macierzy B , natomiast kolumna i pobrana jest z macierzy A. 


[ae oe x x, 
O O a xO u x 
Cae se wt Sk oe (4.14) 
i ee ee a, u 
[a D a O n M 


3. Wyznaczenie wektorów odpowiedzi modelu dla danych zawartych w macierzach: 


A , Bi C. Każdy z wektorów posiada wymiar N x1. 
¥,= f(A), ys=f(B), yc=/(C) (4.15) 
4. Oszacowanie indeksów wrażliwości pierwszego rzędu S, oraz indeksów 
wrażliwości całkowitej S„. Indeks wrażliwości pierwszego rzędu (S,) wyraża 


bezpośredni wpływ parametru i na odpowiedż modelu (indeksy wrażliwości 


wyższych rzędów wyrażają wpływ interakcji między parametrami na odpowiedż 
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układu), natomiast indeks wrażliwości całkowitej (S,, ) określa całkowity wpływ 


parametru i oraz jego interakcji z innymi parametrami na odpowiedź modelu. 


Wskaźniki wrażliwości: pierwszego rzędu oraz całkowitej wyznaczyć można 


korzystając z następujących formuł: 


. yP y% 
FEG] yeng UZRA- 


| P = (4.16) 
V(Y) Va Va-Sto EOT 
‘ 1/N c yy o j 
JIECIA.)] | jadą do 4 Pa = (4.17) 


V(Y) eu, WTZ. 


gdzie: 


s- $e) (4.18) 


Interpretacja uzyskanych wyników. Podczas interpretacji uzyskanych wartości 
indeksów wrażliwości należy zwrócić uwagę na ich następujące właściwości 


[Saltelli i inni, 2004]: 


= wartość indeksu wrażliwości całkowitej S,, powinna być większa niż wartość 


indeksu wrażliwości pierwszego rzędu S,. Możliwa jest sytuacja, w której 
indeksy te są sobie równe - dzieje się tak w jeżeli nie istnieje żadna interakcja 
pomiędzy zmienną X, i pozostałymi czynnikami. Różnica S,,—S, jest miarą 
tego, jak bardzo zmienna X, wplątana jest w interakcję z którąkolwiek z 
pozostałych zmiennych; 

wartość indeksu S,, =0 oznacza, że zmienna X, nie ma wpływu wariancję 
odpowiedzi układu niezależnie od tego, z którymi zmiennymi wchodzi w 
interakcję; 

suma wszystkich indeksów S, jest równa | w przypadku modeli addytywnych 
lub mniejsza od 1 w przypadku modeli nieaddytywnych (zgodnie z definicją z 
modelem addytywnym mamy do czynienia wtedy, kiedy w ramach 


dekompozycji wariancji możliwe jest oddzielenie efektu jego zmiennych 
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wejściowych, np. model Y, 22. jest przykładem modelu addytywnego, 
natomiast Y, = [LZ ; Jest przykładem modelu nieaddytywnego); 


= suma wszystkich indeksów wrażliwości całkowitej jest zawsze większa niż 1. 


Wyjątkowo może być równa 1 w przypadku modeli idealnie addytywnych. 
Do zalet metody badania wrażliwości na podstawie wskaźników Sobola można zaliczyć: 


=" brak zależności między uzyskanymi wynikami a poziomem skomplikowania 
modelu, 

= możliwość uchwycenia wpływu pełnego zakresu zmiany każdego z czynników 
wejściowych na odpowiedź układu, 

= możliwość uwzględnienia w obliczeniach wpływu interakcji pomiędzy 
poszczególnymi czynnikami, 

= możliwość łączenia czynników w grupy i badania wariancji (niepewności jaką dana 


grupa czynników wywołuje) grup czynników. 


Do wad metod zaliczyć należy koszt obliczeniowy - w metodzie MC liczba obliczeń 
potrzebna do wyznaczenia indeksów wrażliwości S, oraz S,,, w przypadku k 
parametrów, wynosi N(k+2). Jeżeli liczba k jest duża (parametrów wejścia jest dużo) 


zaleca się przed właściwą analizą zastosować techniki przesiewowe (Screening Designs). 


N/N MOST 


Rozdział 5. Probabilistyczna analiza wrażliwości - przykłady 


numeryczne 


W niniejszym rozdziale przedstawiono kilka przykładów wykorzystania 
omówionych wcześniej metod i algorytmów. W każdym z nich starano się zaprezentować 
różne aspekty analizy probabilistycznej, uzupełnionej o szacowanie wrażliwości 
spowodowanej zmianą opisujących konstrukcję losowych parametrów. Przykłady 
zaprezentowano w określonej kolejności, począwszy od prostych systemów, a kończąc 
na modelach rzeczywistych konstrukcji. Pomimo tego, iż prezentowane metody mają 
charakter uniwersalny, w pracy ograniczono się jedynie do analizy systemów kratowych. 
Ograniczenie to wynika wyłącznie z nakładu czasu obliczeniowego. W przypadku 
testowania algorytmów probabilistycznych, konieczne jest porównanie kilku metod, a 
więc wykonanie obliczeń dla kilkudziesięciu lub nawet kilkuset realizacji wektora 
zmiennych losowych. Tego typu rozbudowane analizy są praktycznie możliwe jedynie 
dla uproszczonych modeli jakimi są kratownice. Należy jednak podkreślić, że w zakresie 
metodologicznym, obliczeniach innych konstrukcji, np. układów cienkościennych, 


powłok czy płyt można wykonać wykorzystując te same algorytmy. 
5.1 Kratownica von Misesa 
5.1.1 Przedstawienie modelu 


Pierwszy przykład dotyczy symetrycznej kratownicy von Misesa, przedstawionej 
na rys. 5.1. Rozpatrywany model został szczegółowo opisany w [Kleiber i inni, 1997]. 


P 
(1) 


p 1/2 > /2 | 


Rysunek 5.1 Model kratownicy von Misesa 
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Punkty podparcia elementów oddalone są od siebie o / = 93,969 m , natomiast obciążony 
węzeł kratownicy znajduje się na wysokości h= 34,202 m. Pole poprzeczne prętów 
wynosi A=1,0 m” a moduł Young’a E =50,0 kN/m”. Kratownica obciążona jest siłą 
pionową przyłożoną w węźle nr 1. Usytuowanie prętów kratownicy opisane zostało za 
pomocą jednego stopnia swobody- kąta obrotu @ [rad] nachylenia pręta względem 
kierunku poziomego. W modelu założono wstępną imperfekcję związaną z przyjęciem 
wysokości punktu styku prętów. Imperfekcja wyrażona została poprzez kąt a, [rad], 
zawarty pomiędzy kierunkami: założonym i rzeczywistym nachylenia prętów do 


poziomu. 


51.2 Analiza wrażliwości kratownicy na podstawie histogramów stanu 


granicznego 


Na podstawie równania równowagi sił w węźle I wyznaczona została zależność 
obciążenia od kąta obrotu elementów: 
P(9) =-2:N:sin(a — (0, + 9)) (5.1) 
W stanie granicznym powyższa zależność przyjmuje postać: 
P(Q,)=-2:E-A-€,-sin(@a—-(@,+9,)) (5.2) 
gdzie p, jest kątem wychylenia w przypadku maksymalnego skrócenia elementów AŻ,,, 
wyrażonego następująco: 
AL, =€,'L (5.3) 
Do dalszych rozważań kąt wstępnej imperfekcji a, , obciążenie P oraz odkształcenie 


pręta w stanie granicznym €, założono jako zmienne losowe. Dla każdej z wymienionych 


zmiennych wygenerowano zbiór N =1000 realizacji zgodnie z rozkładem normalnym. 


Wartości średnie oraz odchylenia standardowe zmiennych P, £, oraz @ wynoszą: 


m, =0,4, Op =0,05 [kN] (5.4) 
me, =0,02, o, = 0,005 [-] (5.5) 
Mg, = 9,00, Gy, =0,00875 [rad] (5.6) 


Procedurę przedstawionej poniżej analizy wrażliwości w uproszczeniu można opisać w 


następujący sposób: 
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Etap 1: Wybór analizowanego parametru: P, €,, a, 
Etap 2: Generacja zbioru N wartości zmiennych zgodnie z rozkładem normalnym. 
Etap 3: Obliczenie dla każdej realizacji zmiennych wartości mnożnika a, [-] 


analizowanego parametru. W zależności od tego, która zmienna jest 


analizowana, mnożnik œ, przyjmuje wartości: 


_ p = ygr _ Ogr 
a=— ,a4,=-—~,a,=—, (5.7) 


dzie: NÍ., €,, oraz o, są wartościami poszczególnych zmiennych w 
gr ygr Ogr 


stanie granicznym, wyznaczonymi dla i-tej realizacji na podstawie formuły 
32 
Etap 4: Analiza wyników: na podstawie wyznaczonego zbioru bezwymiarowych 


mnożników œ, tworzony jest histogram stanu granicznego. W przypadku 
analizy wrażliwości obciążenia P oraz kąta imperfekcji a, pole histogramu 
dla a,<1,0 odpowiada estymatorowi prawdopodobieństwa awarii p,, który 


wyraża formuła: 
A I 
P, =>} Ia) (5.8) 
N 
gdzie: 


1 d .<1,0 
Ha) =| adi (5.9) 


0 gdy a>10 
W przypadku odkształcenia pręta w stanie granicznym £, estymatorowi 
prawdopodobienstwa awarii p, odpowiada pole histogramu dla a,>1,0. 


Wynik przeprowadzonej analizy w postaci histogramu stanu granicznego 
konstrukcji dostarcza probabilistycznej informacji o zachowaniu się układu 


przy uwzględnieniu losowego charakteru jego parametrów projektowych. 
Pierwszą analizowaną zmienną było obciążenie P . Wygenerowano zbiory 1000 
próbek kąta œ oraz odkształcenia £, . Przy założeniu stanu granicznego konstrukcji, dla 
każdej i-tej realizacji zmiennych £, i œ , obliczono wartość obciążenia granicznego M 7 


według formuły 5.2, a następnie po wygenerowaniu zbioru obciążeń losowych P, dla 


U 
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każdej realizacji wyznaczono mnożnik a, zgodnie z 5.7. Histogram stanu granicznego 
konstrukcji przedstawia rysunek 5.2. Niezawodność ( R = 1— P, ) wyznaczone zgodnie z: 


5.8 15.9 wyniosła R=0,913. 


100 5 


N=1000 realizacji 
70 4 Myi=1,4978; Ogj=0,3771; v=25,18% 
60 4 R=0,913 


częstość [-] 


= 
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2,4 


= 
© 
Rysunek 5.2 Histogram stanu granicznego - analiza wrażliwości zmiennej obciążenia P 


Analogiczne obliczenia wykonano dla zmiennych początkowego kąta a, oraz 
odkształcenia w stanie granicznym €,. Rezultaty — histogramy mnożnika a, 


przedstawiono na rys. 5.3 oraz 5.4. 


18 
16 
14 
= N=1000 realizacji 
3 p mu=57,956; 0,;=187,133; v=28,216% 
2 R=0,978 
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Rysunek 5.3 Histogram stanu granicznego - analiza wrażliwości zmiennej kąta imperfekcji a, 
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oo 
© 


N=1000 realizacji 


4 
oO 


Mai=0,723; Oyi=0,204; v=322,8% 
R=0,907 


częstość [-] 
N wo R a Q 
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Rysunek 5.4 Histogram stanu granicznego - analiza wrażliwości zmiennej odkształcenia €, 


Na podstawie wyników histogramów zmiennej œ, w każdym z rozpatrywanych 
przypadków wyznaczono aproksymowaną wielomianem 2 stopnia funkcję w obszarze 
stanu granicznego a następnie wyznaczono moduł kąta nachylenia lgl stycznej do 
wykresu w punkcie œ, =1,0. Wartości te, uszeregowane malejąco, dostarczają informacji 
o wkładzie analizowanych zmiennych w zmianę odpowiedzi układu. 
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Rysunek 5.5 Obciążenie P- kąt nachylenia stycznej do wykresu w punkcie a, =1 - lol =89°46' 
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Rysunek 5.7 Odkształcenie e, - kąt nachylenia stycznej do wykresu w punkcie 0,=1 - 


|| =89°55" 


Dodatkowej ilościowej informacji o wpływie zmiany analizowanego czynnika na 
odpowiedź układu dostarczają wskaźniki wrażliwości: pierwszego rzędu oraz całkowitej 
Sobola wyznaczone na podstawie wzorów: 4.16 oraz 4.17. W analizowanym przypadku, 
po wygenerowaniu miliona realizacji wyniosły one odpowiednio: Sp =0,1866, 
Srp =0,1943, So, = 0,7934, Se, =0,8196, So, =0,0024, Sræo =0,0055. Otrzymane 
wartości wskazują, że największy wpływ na zmianę odpowiedzi konstrukcji ma zmienna 
odkształcenia e, , natomiast wpływ imperfekcji początkowej (kąta wychylenia prętów 
Œ ) jest znikomy. Te same wnioski można sformułować na podstawie histogramów 


przedstawionych na rys. 5.2, 5.3 oraz 5.4. 
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W prezentowanym przykładzie uzyskano pełne rozkłady stanu graniczne 
uzależnione od przyjętych zmiennych opisujących model kratownicy. Na ich podstawie 
wyznaczono wrażliwość modelu na zmianę poszczególnych parametrów. Opisujące te 
stany zmienne mogą mieć charakter globalny (wartość oczekiwana oraz odchylenia 
standardowe rozkładów stanu granicznego) jak i lokalne — dotyczące punktu 
projektowego funkcji stanu granicznego. Dodatkowych informacji dostarczają także 
wskaźniki wrażliwości Sobola. Należy jednak podkreślić, że pełen opis — histogramy 
stanu granicznego są możliwe tylko w przypadku stosunkowo prostych problemów 


obliczeniowych. 
5.2 Kolumna Zieglera 


Celem analizy jest określenie wrażliwości odpowiedzi modelu na zmianę jej 
parametrów projektowych, opisanych zmiennymi losowymi. Kolumna analizowana była 
w wielu pracach, np. [Alibrandi i inni, 2009; Challamel i inni 2010; Winkelmann, 2013]. 
Przy określaniu wrażliwości kolumny na zmianę jej parametrów projektowych przyjęto 
trzy ścieżki postępowania: 

1. Określenie wrażliwości modelu na podstawie tzw. chmur rozproszenia wyników 
oraz porównaniu współczynników korelacji pomiędzy analizowaną zmienną a 
odpowiedzią układu. 

2. Wyznaczenie wskaźników wrażliwości: pierwszego rzędu oraz całkowitej Sobola. 

3. Określenie wskaźników wrażliwości na podstawie krzywych będących przecięciem 
powierzchni odpowiedzi układu z płaszczyznami zdefiniowanymi poprzez 
przyjęcie wartości badanych zmiennych jako stałe. Powierzchnia mechanicznej 
odpowiedzi modelu wyznaczona została na podstawie punktów pochodzących z 
próbkowania ukierunkowanego TRS opisanego w rozdziale 3.4.3 niniejszej pracy. 
Na potrzebę analizy stworzona została procedura w środowisku MATLAB 
umożliwiająca pozyskanie próbek znajdujących się w bezpośrednim sąsiedztwie 


stanu granicznego g(k,,k,)=0. 


5.2.1 Przedstawienie modelu 


Analizowanym modelem jest utwierdzony u podstawy słup obciążony osiową siłą 


AP. W konstrukcji znajdują się dwa przeguby wewnętrzne (u podstawy i w środku 
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długości słupa) wzmocnione sprężynami o stałych sprężystości k, i k, (rys. 5.8). 
Przyjęto, że obydwie wielkości są zmiennymi losowymi: k, =k,(1+@,), k, =k,(l+a,) 
Zarówno a, jak i a, opisano rozkładem normalnym o wartości średniej równej 


m, =0,0 i odchyleniu standardowym równym Oy, 02. 


Rysunek 5.8 Model kolumny Zieglera [Alibrandi i inni, 2009] 


Funkcja graniczna ma w tym przypadku postać jawną i opisana jest następującym 


równaniem: 


zk (1+ 0)+% (1+a,)-— E (1+) +4kż(1+a,) -A (5.10) 


ej (a, a, ) = 
Przy założeniu, że k,=k,=k,=k,=1[kNm/rad], L=lm oraz P=1kN mnożnik 
obciążenia powodującego wyboczenie konstrukcji wyniósł A=4,=0.5 (3-45 ) ; 


Ostatecznie przyjęto, że A = 0,754; = 0,375 (3 —J5 | oraz założono, że zmienne o i a, 


nie są skorelowane. 
5.2.2 Określenie wrażliwości kolumny na podstawie analizy rozproszenia wyników 


Pierwsza z metod, umożliwiająca szybką i łatwą w przeprowadzeniu analizę 
zachowania modelu jest technika tzw. obrazów rozproszenia, opisana w rozdziale 4.3.1 
niniejszej pracy. Na rysunku 5.9 przedstawiono zależności pomiędzy zmienną 
odpowiedzi układu g a sztywnościami: k, i k,. W przypadku każdej z chmur punktów 


wyznaczono estymator współczynnika korelacji liniowej opisany zależnością: 


p os ,m=1,2 (5.11) 
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Wartości współczynników w przypadku chmur rozproszenia zmiennych k, i k, wyniosły 


odpowiednio: Fxkig=0,899 oraz Fk2g=0,393. 
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Rysunek 5.9 Kształt chmury rozproszenia wyników wyznaczone dla zmiennych k, i k, 


Na rysunku 5.9 czerwonym kolorem zaznaczono wartości Średnie zmiennej 
odpowiedzi układu wyznaczone na podstawie próbek pochodzących z wydzielonych 
zakresów analizowanej sztywności. Zakres zmienności sztywności podzielono na 8 
równych przedziałów. Nachylenie poprowadzonej na podstawie tych punktów prostej 
aproksymacji może stanowić wizualną metodę określenia wrażliwości modelu — wraz ze 
wzrostem kąta nachylenia prostej wzrasta wrażliwość układu na zmianę rozpatrywanego 
parametru. 

Na podstawie wyznaczonych chmur oraz współczynników korelacji można 
wnioskować, że zmienną, na modyfikację której układ jest bardziej wrażliwy jest 


sztywność k, - chmura punktów na pierwszym wykresie „charakteryzuje się stosunkowo 
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łatwym w interpretacji kształtem”, podczas gdy rozkład punktów na drugim wykresie jest 


bardziej zbliżony do jednorodnego. 


5.2.3 Określenie wrażliwości kolumny na podstawie wskaźników wrażliwości 


Sobola 


Kolejnym krokiem analizy jest określenie wskaźników (indeksów) wrażliwości 


Sobola: pierwszego rzędu (S,) oraz wrażliwości całkowitej (Sz) wyznaczone dla 


rozpatrywanych sztywności k, i k,. Indeks wrażliwości pierwszego rzędu S, wyraża 
bezpośredni wpływ i- tego parametru na odpowiedź modelu, natomiast indeks 
wrażliwości całkowitej S$, określa całkowity wpływ i- tego parametru oraz jego 
interakcji z innymi parametrami na odpowiedź końcową. Im większa wartość indeksu 
wrażliwości tym wpływ rozpatrywanego parametru na odpowiedź układu jest większy. 


Wymienione wskaźniki wrażliwości zmiennych k, i k,, wyznaczone na podstawie 


procedur 4.16 - 4.17, wyniosły odpowiednio: S} = 0,7974, ST = (0,8369, S, = 0,1898 


ST, 


az 0,2292. Uzyskane wyniki świadczą o znacząco większym wpływie sztywności 
sprężyny przypodporowej k, na odpowiedź mechaniczną układu, a więc zostały 


potwierdzone wnioski sformułowane na podstawie analizy chmur rozproszenia. 


5.2.4 Określenie wrażliwości kolumny na podstawie analizy krawędzi przecięcia 


powierzchni odpowiedzi z płaszczyznami kig i k2g 


Trzeci sposób określenia wrażliwości zakłada wyznaczenie powierzchni 


odpowiedzi układu $(k) a następnie wyznaczeniu krzywych przecięcia tej powierzchni 
z płaszczyznami k,g oraz k,g. Powierzchnię odpowiedzi aproksymowano na podstawie 


punktów pochodzących z próbkowania ukierunkowanego, a więc takich, które znajdują 
się w bezpośrednim sąsiedztwie stanu granicznego. Na potrzeby analizy stworzono w 
środowisku MATLAB procedurę, która po wpisaniu zakresu oraz liczby kolejnych 
podziałów przestrzeni zmiennych znajduje współrzędne punktów próbkowania. Należy 
podkreślić, że przedstawiona analiza ma charakter jedynie jakościowy — pozwala 


uszeregować zmienne zadania według ich wpływu na odpowiedź układu. 
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Pierwszym etapem analizy było sprawdzenie efektywności próbkowania 
ukierunkowanego. W tym celu porównano ze sobą rezultaty (wartości wskaźnika 
niezawodności Hasofera — Linda /5,, wyznaczone na podstawie wzoru 2.12) otrzymane 
na drodze próbkowania TRS z wynikami uzyskanymi poprzez zastosowanie technik: 
próbkowania: warstwowego (MC-SS) oraz hipersześcianu łacińskiego (MC-LHS). 

Wskaźniki niezawodności /8,, wyznaczono zarówno stosując metodę Monte 
Carlo poszerzoną o w/w techniki redukcji populacji jak również wykorzystując 


powierzchnie odpowiedzi zbudowane na podstawie wyników tych symulacji. 
Metoda Monte Carlo 


Prawdopodobieństwo awarii oszacowane bezpośrednią metodą symulacji Monte 
Carlo, na podstawie 10° probek, wyniosto P; =0,0737. Wskaznik niezawodności, 
wyznaczony zgodnie z formułą 2.8, wyniósł Ø =1,449. 

W przypadku próbkowania warstwowego Monte Carlo obliczenia wykonano dla 
podziału przestrzeni zmiennych na : 2”, 7, e, 5,6, 7,8,9, 10,15, 20°, 25° 
oraz 30° warstw o jednakowych wymiarach. Wartości wskaźników niezawodności 


wyznaczone dla założonych wariantów podziału przestrzeni zawarto w tabeli 5.1. 


W próbkowaniu Monte Carlo poszerzonym o technikę hipersześcianu łacińskiego 
przestrzeń zmiennych została podzielona na 2°, 3°, 4, 5°, 6, 7, 8°, 9°, 10”, 15°, 
20° oraz 25° warstw o jednakowych wymiarach. Losowanie wartości zmiennych 
przebiegało w ten sposób, iż z każdego pojedynczego pasma osi zmiennej k, i każdego 
pojedynczego pasma osi k, dobrana mogła być tylko jedna próbka. Konfiguracja 
rozkładu pobranych w ten sposób punktów utworzona jest w sposób losowy. Wartości 
wskaźników niezawodności Hasofera - Linda uzyskanych dla poszczególnych 


przypadków zagęszczenia warstw w przestrzeni zmiennych zamieszczono w tabeli 5.1. 
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Technika probkowania 
Monte Carlo + SS Monte Carlo + LHS 


Nsauet B [-] 


es 6 [1538 
16 7 


Coa [- | 
5 15508 | 8 | 10 | 
36158001 | o 12445 | 
[6 | 166600 o | - |] 
[9001295158 | 100 | 125651 | 
= 0 __[ 12556 | 
o o f oo 128551 | 
e ee 
[50 E 
= [ET w | 12% | 


Tablica 5.1 Wartości wskaźników niezawodności Hasofera-Linda wyznaczone na podstawie 
bezpośredniego próbkowania Monte Carlo (MC), metody Monte Carlo rozszerzonej o techniki 
próbkowania: warstwowego (MC - SS) i hipersześcianem łacińskim (MC - LHS) 


Metoda Powierzchni Odpowiedzi 


Przy założeniu, że krzywizna rzeczywistej powierzchni odpowiedzi konstrukcji 
jest duża, w obliczeniach wykorzystuje się model drugiego rzędu z członami 
interakcyjnymi, opisanego w przypadku dwóch zmiennych następująco: 

$(k)=B +Bk, + B,k, +B k? + Bk; + By kk, (5.12) 

Aproksymowaną powierzchnię odpowiedzi wyznaczano na podstawie punktów 
obliczeniowych uzyskanych na drodze próbkowania Monte Carlo (zarówno 
bezpośredniego jak również poszerzonego o techniki MC- SS oraz MC- LHS). W celu 
wyznaczenia powierzchni odpowiedzi użyto programu RSM-Win, którego algorytm 
został szerzej opisany w pracy [Winkelmann, 2013]. 

W przypadku bezpośredniego próbkowania Monte Carlo powierzchnia odpowiedzi 
została wyznaczona w 19 wariantach , przy użyciu 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 
200, 300, 400, 500, 600, 700, 800, 900 i 1000 punktów obliczeniowych. Wartości 
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współczynników wielomianu 5.12 oraz wskaźników niezawodności Hasofera- Linda 


wyznaczonych dla każdego przypadku zestawiono w tabeli 5.2. 


10 -0,277 0,272 0,086 -0,091 -0,080 0,184 1,54489 
20 -0,295 0,311 0,084 -0,096 -0,064 0,155 1,53122 
30 -0,354 0,294 0,211 -0,090 -0,125 0,161 1,54149 
40 -0,344 0,305 0,175 -0,091 -0,101 0,151 1,52312 
50 -0,343 0,306 0,174 -0,091 -0,100 0,150 1,52484 
60 -0,343 0,305 0,176 -0,090 -0,101 0,149 1,52610 
70 -0,335 0,303 0,165 -0,092 -0,099 0,153 1,53774 
80 -0,317 0,279 0,152 -0,091 -0,102 0,175 1,54716 
90 -0,317 0,279 0,152 -0,091 -0,102 0,174 1,54712 
100 -0,319 0,281 0,155 -0,091 -0,103 0,173 1,54830 
200 -0,301 0,270 0,134 -0,089 -0,098 0,180 1,55492 


1,55249 
1,55375 
1,55047 
1,54962 


1,54602 
1,54704 
1,54684 
1,54656 


Tablica 5.2 Współczynniki powierzchni odpowiedzi wyznaczone na podstawie zbioru N «4mp: 


punktów obliczeniowych (próbkowanie Monte Carlo) oraz wyznaczone na ich podstawie wskaźniki 
niezawodności Hasofera-Linda 


Aproksymacja powierzchni w przypadku punktów obliczeniowych pochodzących 
z próbkowania warstwowego Monte Carlo wykonana została dla 12 wariantów podziału 
przestrzeni zmiennych. Wartości współczynników regresji dla poszczególnych 


przypadków podziału oraz uzyskane wskaźniki niezawodności /8,, zamieszczono w 


tabeli 5.3. 

Podobne obliczenia wykonano w przypadku powierzchni odpowiedzi 
wyznaczonej na podstawie punktów uzyskanych z próbkowania Monte Carlo 
poszerzonego o technikę hipersześcianu łacińskiego. W tym wypadku przestrzeń 
zmiennych została podzielona na 2.3.4.5,6 .,7.,8.,9.10 157,20 oraz 25° 
warstw o jednakowych wymiarach. Zestawienie wyników dla tego typu próbkowania 


zawiera tabela 5.4. 
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| 9 | -0,287 0,271 0,111 | -0,089 -0,091 0,180 | 1,55537 | 


16 -0,289 0,271 0,114 -0,089 -0,092 0,181 1,55559 
25 -0,290 0,270 0,115 -0,089 -0,093 0,181 1,55570 
36 -0,290 0,270 0,117 -0,089 -0,093 0,181 1,55567 


81 -0,292 0,269 0,120 -0,089 -0,095 0,181 1,55545 
100 -0,293 0,269 0,121 -0,089 -0,095 0,182 1,55542 
225 -0,294 0,269 0,124 -0,089 -0,096 0,182 1,55505 
400 -0,295 0,269 0,125 -0,089 -0,097 0,182 1,55484 
625 -0,296 0,269 0,127 -0,089 -0,097 0,182 1,55467 


1,55455 


Tablica 5.3 Współczynniki powierzchni odpowiedzi wyznaczone na podstawie zbioru N 4mp: 


punktów obliczeniowych pochodzących z próbkowania warstwowego Monte Carlo oraz 
wyznaczone na ich podstawie wskaźniki niezawodności Hasofera-Linda 


Na rys 5.10 przedstawiono graficzne porównanie uzyskanych zbieżności 
iteracyjnych obliczeń wskaźnika niezawodności /,, w przypadku powierzchni 
odpowiedzi aproksymowanej na podstawie punktów pochodzących z wyżej opisanych 
rodzajów próbkowania. Uzyskane wyniki wskazują, że w przypadku analizowanego 
modelu za pomocą metody Monte Carlo, przy 1000 realizacji, nie otrzymano 
ustabilizowanego rozwiązania. Spośród analiz wykorzystujących metodę RSM 
próbkowanie SS pozwala na uzyskanie wiarygodnych wyników już na poziomie kilku 
próbek. We wszystkich pozostałych przypadkach analizy RSM, zwiększanie liczby 
próbek nie prowadzi to stabilizacji rozwiązania. Należy jednak podkreślić, że różnice 
wyników dotyczą drugiego lub nawet trzeciego miejsca po przecinku, a więc zastosowane 
różne podejścia prowadzą do rozwiązań, które można generalnie określić jako 


powtarzalne i ustabilizowane z inżynierskiego punktu widzenia. 
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0,00000 
1,65922 
0,87465 
1,55090 
1,56646 
1,54398 
1,56939 
1,55393 
1,54111 
1,54945 
1,55845 
1,54712 
1,54281 
1,54701 
1,55280 
1,55557 
1,55170 


1,55087 
1,55767 


Tablica 5.4 Współczynniki powierzchni odpowiedzi wyznaczone na podstawie zbioru N «4mp: 


punktów obliczeniowych (próbkowanie Monte Carlo wraz z techniką hipersześcianu łacińskiego) 
oraz wyznaczone na ich podstawie wskaźniki niezawodności Hasofera-Linda 
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Rysunek 5.10 Zbieżność obliczeń wskaźnika niezawodności /,, oszacowanego w przypadku 


bezpośredniego próbkowaniu Monte Carlo a także na podstawie powierzchni odpowiedzi 
aproksymowanych przy użyciu punktów pochodzących z 4 różnych technik próbkowania 


N/N MOST 


N/N MOST 


5.2 Kolumna Zieglera 17 


Technika próbkowania ukierunkowanego (TRS) 


W celu wyznaczenia punktów próbkowania TRS, w środowisku MATLAB 
stworzono procedurę, która w sposób automatyczny zagęszcza punkty w przestrzeni 
zmiennych a następnie z całego zbioru wybiera te, które znajdują się w pobliżu stanu 
granicznego. Struktura procedury zbudowana została na podstawie algorytmu 


przedstawionego w rozdziale 3.4.3 pracy. Całość składa się z następujących modułów: 


1. Wczytanie danych wejściowych: liczby zmiennych zadania, maksymalnej liczby 
podziałów przestrzeni podczas której poszukiwana jest para warstw, z której jedna 
zawiera zdarzenie f (takie dla którego funkcja stanu granicznego g(x) <0 ), druga 
zaś zdarzenie s (takie dla którego zachodzi g(x) > 0 ), ostatecznej liczby podziałów 
przestrzeni Q , wektora zawierającego zakresy dla poszczególnych sztywności oraz 


funkcji stanu granicznego. Przyjęto, że sztywności k, i k, mogą przyjmować 
wartości z zakresu (k; —30;k, +30) a zatem (0,4;1,2) [kNm/m]. 


2. Pierwotny podział przestrzeni przy użyciu techniki RSS, do momentu uzyskania 
minimum jednej pary sąsiadujących ze sobą warstw zawierających zdarzenia f - s. 
W przypadku gdy program osiągnie liczbę podziałów przyjętą przez użytkownika 
jako maksymalną, nie znalazłszy pary f-s wyświetlony zostaje stosowny 
komunikat o konieczności zwiększenia liczby podziałów. Losowania próbki z 
wydzielonej warstwy odbywa się zgodnie z rozkładem równomiernym. 

3. Wybór spośród wyłonionych par f-s tej, której waga jest największa 
Wax = Wy + W, Przez wagę danej warstwy (w, lub w, ) określa się iloczyn długości 
jej boków we wszystkich wymiarach (np. w problemie dwuwymiarowym będzie to 
pole danej warstwy). Jeżeli więcej niż jedna para warstw ma taką samą wagę w,,,. 
następuje losowy wybór jednej z par. 

4. Wyznaczenie długości odcinka łączącego próbki z danej pary warstw oraz jego 
dłuższego rzutu (zgodnie z formułami: 3.39 oraz 3.40). W przypadku wystąpienia 


kilku rzutów o tej samej długości wybierany jest losowo jeden z nich. 


5. Wyznaczenie współrzędnej |x 


cut 


| (rys.3.5) - punktu przecięcia prostej prostopadłej 


do dłuższego rzutu odcinka łączącego próbki w warstwach /-s, zgodnie z formułą 


(3.41). 
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6. Sprawdzenie w której spośród warstw (f czy s) nastąpi podział. W przypadku gdy 
punkt przecięcia pokrywa się z krawędzią podziału warstw program jako miejsce 
podziału przyjmuje połowę dłuższego rzutu odcinka łączącego próbki f— s. 

7. Kolejno wykonywane są: podział warstwy, wylosowanie nowej próbki oraz 
powtórzenie całej procedury zadaną liczbę razy. 

8. Wyznaczenie prawdopodobieństwa awarii rozumianego jako objętość warstw 
zawierających zdarzenie awaryjne (wzór 3.42). 

9. Utworzenie pliku zawierającego współrzędne uzyskanych próbek oraz wartość 


funkcji stanu granicznego wyznaczoną dla tych współrzędnych. 


Wykorzystując procedurę TRS przestrzeń zmiennych podzielono na warstwy, z 


których wylosowano po jednej próbce. 


Przykłady wygenerowanych punktów 


obliczeniowych w przypadku pośrednich podziałów na 50 i 200 warstw oraz ostatecznego 


podziału na 1000 warstw przedstawiono na rys. 5.11. 
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Rysunek 5.11 Punkty obliczeniowe pochodzace z: 50, 200 oraz 1000 warstw wyznaczonych na 
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podstawie procedury TRS 
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Na podstawie współrzędnych uzyskanych na drodze kolejnych podziałów 
aproksymowano powierzchnię odpowiedzi $ (k,,k, ). Wartości współczynników regresji 


dla kolejnych podziałów przestrzeni oraz uzyskane wskaźniki niezawodności Hasofera — 


Linda zamieszczono w tabeli 5.5. 


Bu 
0,250 | -0,750 | 0,625 1,688 | -0,219 | -0,555 
1,73232 
-0,40183 
1,59534 
1,59281 
1,57906 
1,58618 
1,57873 


1,56169 
1,56320 
1,56720 
1,57278 
1,56289 
1,56462 
1,56159 
1,55890 
1,56100 
1,56129 
1,56223 
1,56494 
1,56340 
1,56172 
1,55971 
1,55913 


Tablica 5.5 Współczynniki powierzchni odpowiedzi wyznaczone na podstawie zbioru N 4mp: 


punktów obliczeniowych (próbkowanie TRS) oraz onliczone na ich podstawie wskaźniki 
niezawodności Hasofera-Linda 


W przypadku analizowanego modelu próbkowanie TRS nie poprawiło 
rozwiązania, szczególnie jeżeli porównać je do wyników uzyskanych za pomocą 


próbkowania warstwowego. Trudno też dostrzec stabilizację wyników w przypadku 


wzrastającej liczby próbek. Jednakże, podobnie jak miało to miejsce we wcześniejszych 
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rozwiązaniach, różnice we wszystkich uzyskanych wynikach są niewielkie i dotyczą 
drugiej lub nawet trzeciej liczby po przecinku. Metodę TRS należy ocenić jako bardziej 
wiarygodną, gdyż generacja punktów dotyczy jedynie obszaru funkcji stanu granicznego 
(rys. 5.11). Należy spodziewać się, iż tego typu analiza będzie miała większe znaczenie 
w przypadku bardziej skomplikowanych modeli, które charakteryzują się znacznie 


większymi rozrzutami rozwiązań. 


Podstawą badanie wrażliwości odpowiedzi układu w tym przypadku jest 
obserwacja przebiegu powierzchni odpowiedzi w poszczególnych obszarach przestrzeni 
zmienności parametrów k, i k,. Szacowania wrażliwości układu dokonano na podstawie 
1000 punktów pochodzących z próbkowania TRS. Równanie powierzchni opisane jest 
następująco: 

&(k,,k,) =-0,166 + 0,162k, -0,048k, - 0,098k; -0,044k; +0,248k,k, (5.13) 
Na rys. 5.13 przedstawiono krawędzie przecięcia powierzchni odpowiedzi z 


płaszczyznami poprowadzonymi dla 4 różnych wartości zmiennych k, i k,. Schemat 


przyjętych płaszczyzn przecięcia zaprezentowano na rys. 5.12. 
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Rysunek 5.12 Schemat przyjętych płaszczyzn przecięcia (linie czerwone, przerywane) powierzchni 
odpowiedzi kolumny 
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Rysunek 5.13 Krawędzie przecięcia powierzchni odpowiedzi układu z płaszczyznami k,g i k,g 


Na podstawie zamieszczonych na rys.5.13 wykresów zaobserwowano, że przy 
stałej wartości zmiennej k, powierzchnia odpowiedzi układu charakteryzuje się 
mniejszym nachyleniem do poziomu. Należy zatem uznać, że jest to parametr, na który 
zmienna wyjściowa jest wrażliwsza. Wnioski te znajdują potwierdzenie we 


wcześniejszych obliczeniach wskaźników wrażliwości Sobola. 
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5.3 Wieża telekomunikacyjna 


W niniejszym rozdziale rozpatrywana jest przestrzenna wieża kratownicowa 
będąca elementem infrastruktury sieci telekomunikacyjnej, opisana szerzej w [Szafran i 


inni, 2016], przedstawiona na rysunku 5.14. 
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Rysunek 5.14 Schemat wieży telekomunikacyjnej - widok z boku (a) oraz z góry (b) (na podstawie 
[Szafran i inni, 2016]) 
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5.3.1 Przedstawienie modelu 


Wieża o całkowitej wysokości 40m, podzielona jest na 7 segmentów. Dolna część 
konstrukcji (segmenty 2-7, znajdujące się poniżej 34m) ma kształt Ściętej piramidy o 
stałym nachyleniu równym 5% . Długość boku podstawy piramidy zmienia się od 4,9m 
do 1,5m. Elementy segmentu nr 1 tworzą pryzmat o długości boku podstawy równej 
1,5 m. Widok wieży przedstawiono na rysunku 5.14. 

Elementy nośne konstrukcji (krawężniki) zostały wykonane z okrągłych prętów 
pełnych, natomiast elementy wykratowania wykonano z gorącowalcowanych 
kątowników równoramiennych i nierównoramiennych. Szczegóły dotyczące przekrojów 
poszczególnych elementów wieży zestawiono w tablicy 5.6 . 

System ukośnych stężeń wieży stanowi wykratowanie typu X. Elementy na 
przecięciu połączone są ze sobą za pomocą odstępnika i pojedynczej śruby, natomiast 
elementy wykratowania połączone są z nogami wieży za pomocą blachy węzłowej oraz 
dwóch śrub. W modelu matematycznym system odstępników i blach węzłowych nie jest 
ujęty — pręty łączą się w węzłach bezpośrednio, tworząc przestrzenne połączenia 
przegubowe. Zestawienie właściwości mechanicznych stali, z której wykonana jest 
wieża, przedstawia tablica 5.7. Podstawą do ich ustalenia były wyniki badań rozciągania 
18 próbek (12 o przekroju kołowym oraz 6 o przekroju kątowym), zaprezentowane w 


[Szafran i inni, 2016]. 


S-1 © 65 L 60x60x5 
S-2 0 65 L 60x60x5 


S-3 Ø 80 L 60x60x6 
S-5 © 90 L 90x60x8, L 100x75x8 
S-6 © 90 L 100x75x8 


S-7 © 100 L 120x80x8 


Tablica 5.6 Profile wybranych elementów wieży (wszystkie wymiary podane są w mm) [Szafran i 
inni, 2016] 


moduł Young’a 202,0 GPa 206,6 GPa 
granica plastyczności 294,2 MPa 273,2 MPa 


wytrzymałość na rozciąganie 399,8 MPa 284,1 MPa 


Tablica 5.7 Właściwości mechaniczne stali konstrukcyjnej elementów składowych wieży [Szafran i 
inni, 2016] 
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Wieża obciążona jest jedną siłą skupioną, modelującą: 
1) wypadkową działania wiatru na wieżę, w przypadku, gdy kąt natarcia oraz wartość 
tego oddziaływania są najmniej korzystne, 


2) wypadkową działania wiatru na elementy wyposażenia teletechnicznego wieży. 


Punkt przyłożenia siły został ustalony po zredukowaniu przedstawionego na 
rysunku 5.15 układu sił do jednej wypadkowej. W pracy [Szafran i inni, 2016] 
udowodniono, iż rezultaty uzyskane na podstawie modelu numerycznego z 
uproszczonym schematem obciążenia są zbieżne z obliczeniami przeprowadzonymi przy 


wykorzystaniu modelu z obciążeniem rozłożonym na wysokości wieży. 


N 


N | 


F=9,77kN : b. 
2x097 Nm | j K SK 

= 11,91 kN pz- í pí 
eona [dg pi 
oom AE l 
A 2) 

L- DI F=116,30 kN IN 

21091 Nim |-- q y | «i ś 
maz. | X 
NK NN 

2x0,90 kN/m | : SK 
W 

Li Ne: 
HA ki 
2x092%Nim | N U NA 
Ba > s GD 
By 4 a AO 
22091 Nim | OH) O| 
XAN my 
l Sh MK 
BT ZY 

Y o a 


Rysunek 5.15 Schemat obciążenia wiatrem wieży telekomunikacyjnej i jej wyposażenia 
teletechnicznego — obciążenie pełne (a) oraz jego redukcja do wypadkowej (b) 
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Wieża posadowiona jest na gruncie w sposób pośredni, za pomocą stalowej ramy 
stelażowej, do której przymocowano węzły podstawy konstrukcji. Rama stalowa 
spoczywa na poduszce z płyt żelbetowych, co ma zapewnić odpowiednie przekazanie 
obciążeń z ramy na podłoże gruntowe. Dodatkowo, rama dociążona jest balastem z płyt 
żelbetowych. 

Badania przeprowadzone w pracy [Szafran i inni, 2016] wykazały jednoznacznie, 
iż sztywność układu: rama — podłoże gruntowe jest istotnym czynnikiem wpływającym 
na stany graniczne: stateczności i użyteczności wieży. Ustalono, że zmiany w sztywności 
podparcia węzłów podstawy wieży doprowadzić mogą do zmniejszenia mnożników: 
dopuszczalnego lub krytycznego obciążenia wieży, zarówno z uwagi na przekroczenie 
dopuszczalnego przemieszczenia jej górnego węzła (u,,,=H/100=40mm) jak 
również z uwagi na zmianę miejsca w konstrukcji, w którym dochodzi do utraty 
stateczności. 

Istotnym problemem staje się wówczas obliczenie parametrów określających 
zmienność sztywności podparcia. Można tego dokonać badając podatność układu: 
rama — podłoże gruntowe na działanie jednostkowego obciążenia: Ściskającego i 
rozciągającego. W tym celu w programie SOFiSTiK utworzono model numeryczny 
układu podpierającego bazując na wymiarach rzeczywistej konstrukcji, zamieszczonych 
w pracy [Szafran i inni, 2016]. Widok rzeczywistej ramy stelażowej — ramy na poduszce 
z płyt żelbetowych z balastem oraz bez przedstawiono na rys 5.16. Wizualizację modelu 
ramy stelażowej przedstawiono na rysunku 5.17. W modelu tym sztywność podłoża 
gruntowego (przyjęto piasek średni, średnio zagęszczony) oraz sztywność poduszki z płyt 
betonowych zostały przeliczone na zastępcze sztywności elementów sprężystych zgodnie 


z procedurami programu SOFiST1K. 


dodatkowe 


(a) (b) 


Rysunek 5.16 Ilustracja rzeczywistej ramy stelażowej podpierającej wieżę telekomunikacyjną 
— rama na poduszce z płyt żelbetowych (a) i dociążenie balastem (b) [Szafran i inni, 2016] 
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<Ś soristik 


Rysunek 5.17 Model ramy stelażowej (a) oraz dociążenia balastem (b) wykonany w programie 
SOFiSTiK 


W modelu ramy przyjęto następujące założenia o możliwych rozrzutach parametrów 
wejściowych modelu: 
= rozrzut współczynnika sprężystości podłoża gruntowego, związana ze stopniem 
zagęszczenia piasku średniego, może wahać się w granicach od 22,5 do 
67,5 MN/m?. Zgodnie z PN-81/B-03020-7 zmienna ta opisana jest rozkładem 
równomiernym; 
= wartość modułu sprężystości stali E,, może wahać się w granicach od 184,5 do 
225,5 GPa. Zmienną opisano rozkładem normalnym; 


= wartość modułu sprężystości betonu C20/25 - E, może wahać się w granicach od 


21 do 39 GPa. Zmienna ta również opisana jest rozkładem normalnym; 
= zmienność obciążenia balastem, opisana rozkładem równomiernym, może wynosić 
do 25% wartości podstawowej. Związane jest to z możliwą lokalną zmianą 


objętości pasma balastu. 


Dla wyżej zaproponowanych zmiennych, na podstawie wartości Średniej 
naprężenia pod poduszką betonową ramy, obliczono zastępcze współczynniki k w dwóch 
wariantach: od jednostkowego obciążenia ściskającego pod słupem ściskanym wieży 
oraz od pary jednostkowych obciążeń rozciągających pod słupami poddanymi 
rozciąganiu. Obliczone wartości zostały wprowadzone jako stałe sprężystości w 
podporach sprężystych, założonych w węzłach podporowych modelu wieży. Dzięki 
takiemu postępowaniu w modelu numerycznym w prosty sposób uwzględniono ramę 
podporową. Wizualizację modelu wieży punktowo sprężyście podpartej przedstawiono 


na rysunku 5.18. 
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Rysunek 5.18 Model rozważanej wieży telekomunikacyjnej 
punktowo sprężyście podpartej. 


Na podstawie obliczeń statycznych modelu przedstawionego na rysunku 5.18, 


wyznaczono następujące wartości zastępczego współczynnika k podpory sprężystej: 


= w przypadku sprężyny pod słupem ściskanym wartość średnia oraz odchylenie 
standardowe współczynnika k, wyniosły odpowiednio: “4, =240 MN/m, 
0,, =65 MN/m, funkcja gęstości prawdopodobieństwa w przybliżeniu opisana 
jest rozkładem normalnym; 

= w przypadku sprężyny pod słupami rozciąganymi — wartość średnia oraz 


odchylenie standardowe współczynnika k, wyniosły: 4 =360 MN/m, 


N/N MOST 


88 Rozdziat 5 Probabilistyczna analiza wrazliwosci - przyktady numeryczne 


o, =100 MN/m, funkcja gęstości prawdopodobieństwa w przybliżeniu opisana 


jest rozkładem normalnym. 


Można więc założyć, iż rozważany w zadaniu problem jest probabilistycznie 
dwuwymiarowy. Uproszczenie takie umożliwia graficzną prezentację zmiany 


odpowiedzi konstrukcji (wyrażonej poprzez mnożnik obciążenia dopuszczalnego m,,, ) 


na skutek wahania wartości współczynników sprężystości podpór k, i k,. 


5.3.2 Wstępna analiza zmiany odpowiedzi konstrukcji 


Pierwszym krokiem analizy modelu było przeprowadzenie tzw. 
eksperymentowania czynnikowego (opisanego w rozdziale 4.3.1 pracy), mającego na 
celu obserwację mechanicznego zachowania odpowiedzi modelu na skutek zmiany 
parametrów kc i kı. Analiza taka ma charakter jedynie szacunkowy, dlatego liczba użytych 
na jej potrzeby realizacji zmiennych powinna być mała. Jednocześnie istotnym jest 
właściwy dobór punktów w przestrzeni realizacji, dla których badana jest odpowiedź 
konstrukcji - niewłaściwie przyjęte punkty prowadzą do przekłamań w obrazie 
zachowania modelu. Właściwie przeprowadzone badania czynnikowego umożliwia 
rozpoznanie jaki charakter - liniowy czy nieliniowy ma odpowiedź konstrukcji oraz gdzie 
w przestrzeni zmiennych orientacyjnie znajduje się stan graniczny. Właściwie 
wyciągnięte wnioski umożliwiają określenie, które z metod próbkowania potencjalnie 
mogą być właściwe w bardziej szczegółowej analizie a które należy odrzucić. 

W zadaniu rozpatrzono dwa warianty zbioru punktów. Pierwszy zawiera dziewięć 
próbek będących wynikami uzyskanymi na podstawie wartości zmiennych 
odpowiadających ich wartościom średnim oraz wartościom średnim pomniejszonym i 
powiększonym o jedno odchylenie standardowe (u—06;u;u+0 ). Dla każdej próbki 
wykonano obliczenia numeryczne określając mnożnik obciążenia dopuszczalnego ze 
względu na kryteria: SGU i SGN. Wyniki obliczeń przedstawiono w tabeli 5.8. Ilustrację 


graficzną śledzonej zmiany zaprezentowano na rysunku 5.19. 
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x, =k, [MN/m] 
U—0=45 U =240 U+0 =435 


1,562 1,654 1,714 


1,602 1,696 1,758 
u+ =660 1,626 1,722 1,785 


Tablica 5.8 Wartości mnożnika obciążenia dopuszczalnego wieży ze względu na kryteria SGU i 
SGN wyznaczone dla dziewięciu próbek . 


2, , < 
“2 16 + y~ = 1,6-2 
£ 12 T 4 N 460 m 1,2-1,6 

0,8 < m 0,8-1,2 

Ng 360 
KW 
ke[MN/m] 240 N ann 
260 


305 


Rysunek 5.19 Zmiana mnoznika obciążenia dopuszczalnego wieży na skutek zmiany parametrów 
k, i k,- dobór zmiennych ze zbioru (4-0; 4; u +0) 


Drugi wariant zakłada przyjęcie dziewięciu próbek zawierających wielkości: 


średnie oraz ekstremalne współczynników odpowiadające u—30; u; +30. Dla każdej 


próbki ponownie wykonano obliczenia, określając mnożnik obciążenia dopuszczalnego. 


Wyniki obliczeń zestawiono w tabeli 5.9. Ilustrację graficzną śledzonej zmiany 


przedstawia rysunek 5.20. 
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Rysunek 5.20 Zmiana mnożnika obciążenia dopuszczalnego myo, wieży na skutek zmiany 


parametrów ką i k, - dobór zmiennych ze zbioru (4 —30; 4; U +30) 


W przypadku pierwszego wariantu doboru punktów przebieg zmienności 
mnożnika dopuszczalnego wieży nie wskazuje na wyraźną miejscową nieliniowość 
zagadnienia (widoczną przy zastosowaniu drugiego wariantu obliczeń). Aproksymacja 
powierzchni odpowiedzi wielomianem stopnia pierwszego w tym przypadku wydaje się 
być wystarczającą. W drugiej wersji doboru punktów, jak można zaobserwować na 
rysunku 5.20, zmiana mnożnika obciążenia dopuszczalnego ma charakter nieliniowy, 
gdyż w przypadkach, w którym zmienne przyjmują minimalne wartości następuje 
zauważalny spadek mnożnika. Całość estymowanej powierzchni zdaje się być możliwa 
do aproksymacji jedynie wielomianami stopnia większego niż 1. Dodatkowo, na 
podstawie danych zebranych w tablicy 5.9 można stwierdzić, iż wieża w analizowanych 
zakresach zmienności przekracza stan graniczny z uwagi na dwa kryteria. Wynik 
oznaczony symbolem (*) wskazuje na sytuację, w której dla małych wartości 
współczynników k istnieje możliwość wystąpienia nadmiernego przechyłu wieży, 
przekraczającego wartość dopuszczalną ze względu na SGU. Wyniki oznaczone 
symbolem (**) sygnalizują możliwość wystąpienia sytuacji, w której przy wysokich 
wartościach współczynników k następuje lokalna utrata stateczności elementów wieży. 


Obserwacja mnożnika m,,, Z uwagi na kryterium stateczności wskazuje, że w przypadku 
małych wartości k mnożnik m,,,, przyjmuje wartości większe (dochodzące do 1,832), a 


utrata stateczności pojawia się w górnych partiach wieży (związane jest to z nadmiernym 
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przechyłem konstrukcji oraz lokalizacją obciążenia zewnętrznego), natomiast dla 
średnich i dużych wartości k mnożniki krytyczne są mniejsze (wahają się w granicach 
1,811 — 1,823), a utrata stateczności pojawia się wówczas w elementach wieży 


znajdujących się blisko punktów podparcia. 
Analiza zmiany odpowiedzi konstrukcji wykorzystująca metodę Monte Carlo 


Kolejny krok analizy zakładał obliczenie mnożnika m,,, w 200 modelach, w 


których współczynniki k przyjmują wartości wygenerowane zgodnie z metodą Monte 


Carlo, a więc całkowicie dowolnie i niezależnie. Otrzymano następujące ekstremalne 


min max 


wartości mnożnika: m, =Vmin =1,119 oraz my = Ymax =1821, z czego druga z 


wymienionych wartości nie spełnia kryterium stateczności. 

Na rysunkach: 5.21 oraz 5.22 przedstawiono proces zmiany wartości średniej oraz 
odchylenia standardowego mnożnika obciążenia dopuszczalnego wieży, a na rysunku 
5.23 poglądowy histogram wyników mnożnika dopuszczalnego dla 200 realizacji 


zmiennych metodą MC. 


wartość średnia mnożnika Mao [-] 


liczba próbek Nsavpt [-] 


1 21 41 61 81 101 121 141 161 181 
Rysunek 5.21 Proces zmiany wartości średniej mnożnika obciążenia dopuszczalnego wieży dla 200 
próbek wygenerowanych według metody Monte Carlo 


N/N MOST 


92 Rozdziat 5 Probabilistyczna analiza wrazliwosci - przyktady numeryczne 


liczba próbek Nsavpt [-] 


odchylenie standardowe mnoznika Mo, [-] 


0 T T T T T T T T T T 
1 21 41 61 81 101 


121 141 161 181 


Rysunek 5.22 Proces zmiany odchylenia standardowego mnożnika obciążenia dopuszczalnego 
wieży dla 200 próbek wygenerowanych według metody Monte Carlo. 


Hm,,, 1,6683 [- 
Omy,,=0,1053 [H 


częstość [-] 


wo 
© 


1 1,1 1,2 1,3 1,4 1,5 1,6 1,7 1,8 


mnożnik dopuszczalnego obciążenia wieży Mggp [-] 


Rysunek 5.23 Histogram mnożnika obciążenia dopuszczalnego wieży utworzony na podstawie 200 
próbek wygenerowanych według metody Monte Carlo. 


Analiza zmiany odpowiedzi konstrukcji wykorzystująca metodę próbkowania 


ukierunkowanego TRS 


W kolejnym kroku analizy, w związku z zauważalnym spadkiem powierzchni 
odpowiedzi, zdecydowano się prześledzić szczegółowo obszar przekroczenia funkcji 
stanu granicznego konstrukcji. W tym celu wykorzystano stworzoną na podstawie 
algorytmu metody TRS autorską procedurę generacji próbek. 

Procedura ta wymaga podania przez użytkownika równania funkcji stanu 
granicznego, w pobliży której mają być losowane próbki. Wykorzystano 15 realizacji 


wektora zmiennych metodą Monte Carlo, na podstawie których w programie RSM - WIN 
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wyznaczono współczynniki regresji powierzchni odpowiedzi konstrukcji. 
Aproksymowana powierzchnia dana jest równaniem: 
P(x, x, )=—0,41779 + 4,4821- x, +1,6654- x, + (5.14) 
— 6,6876: x; —1,8618- x; —1,059- xx, 


natomiast równanie jej przecięcia z płaszczyzną f(x,,x,) = y(x,,x,)=1 opisuje formuła: 
x, (x, ) = 268557 -(*), (5.15) 


gdzie: 


(*) = 0,00167-1,059-10 x, —6,98-10* -|--6939,03 +613,19x, -x (5.16) 


Przebieg przybliżonej funkcji stanu granicznego oraz rezultat generacji techniką 


próbkowania ukierunkowanego przedstawiono na rysunku 5.24. 


200 4 


— przybliżona funkcja stanu granicznego 
(uzyskana na podstawie 15 próbek MC) 
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Rysunek 5.24 Przybliżona funkcja stanu granicznego wyznaczona na podstawie 15 próbek symulacji 
MC oraz 200 punktów pochodzących z próbkowania techniką TRS 


Podobnie jak w przypadku metody Monte Carlo w procesie TRS wygenerowano 
200 próbek (modeli numerycznych). Otrzymano następujące wartości ekstremalne 


mnożnika obciążenia dopuszczalnego: min =0,992 oraz mj; =1,225, z czego 


pierwsza z wymienionych wartości, jak większość w tej analizie, nie spełnia kryterium 
SGU. Rysunek 5.25 przedstawia poglądowy histogram utworzony na podstawie wyników 
mnożnika obciążenia dopuszczalnego dla 200 realizacji pochodzących z próbkowania 


TRS. 
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Rysunek 5.25 Histogram mnożnika obciążenia dopuszczalnego wieży utworzony na podstawie 200 
próbek wygenerowanych zgodnie z metodologią TRS 


Analiza zmiany odpowiedzi konstrukcji wykorzystująca metodę powierzchni 


odpowiedzi (RSM) 


Uzyskane w poprzednich punktach analizy następujące próbki: 1 próbka 
zawierająca wartości średnie zmiennych, 8 próbek zawierających wielkości ekstremalne 
(odpowiadające wut30 ) zmiennych, 8 próbek odpowiadających “+o, 200 próbek 
symulacji Monte Carlo oraz 200 próbek techniki symulacji TRS zostały w tym kroku 
wykorzystane do określenia współczynników kierunkowych aproksymowanej 
powierzchni odpowiedzi konstrukcji 

Z uwagi na widoczną nieliniowość pojawiającą się w obrębie „ogonów” 
rozkładów prawdopodobieństwa obydwu zmiennych, aproksymacja wykonana zostanie 
wielomianem drugiego stopnia. Równanie powierzchni odpowiedzi opisuje wówczas 
następująca formuła: 

P(x, x) = B, + Bx, + B,x, +B x, + Bx) + ByXX (5.17) 
Należy podkreślić, ze liczba próbek w zbiorach celowych (9 współrzędnych zmiennych 
x, 1 x,) jest zbyt mała aby właściwie aproksymować powierzchnię odpowiedzi 


konstrukcji. Zdecydowano zatem aby próbki te zostały włączone do zbioru wyników 
pochodzących z symulacji Monte Carlo. Ostatecznie przyjęto 4 warianty obliczeniowe - 


ich opis zawarty jest w tablicy 5.10. 
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NR WARIANTU 
OBLICZENIOWEGO 


SPOSOB DOBORU PROBEK 


1 próbka zawierająca wartości Średnie, 8 próbek zawierających 
wartości ekstremalne (odpowiadające u+ 30 ), 200 próbek symulacji 
MC — łącznie 209 próbek; aproksymacja powierzchni odpowiedzi po 
9 pierwszych i każdych 5 kolejnych próbkach; 


1 próbka zawierająca wartości Średnie, 8 próbek zawierających 


wartości odpowiadające “+o, 200 próbek symulacji MC- łącznie 


209 próbek; aproksymacja powierzchni odpowiedzi po 9 pierwszych i 


każdych 5 kolejnych próbkach; 


200 próbek symulacji MC — aproksymacja powierzchni odpowiedzi po 
każdych 5 kolejnych próbkach; 


200 próbek symulacji TRS — aproksymacja powierzchni odpowiedzi 
po każdych 5 kolejnych próbkach. 


Tablica 5.10 Wartości mnożnika obciążenia dopuszczalnego wieży ze względu na kryteria SGU i 
SGN dla dziewięciu próbek dobranych wg metodologii PEM 


Schemat doboru punktów w kolejnych wariantach obliczeniowych przedstawiono 


na rys. 5.26. Na rysunku 5.27 zaprezentowano zestawienie próbek wygenerowanych 


zgodnie z metodami: MC (wariant 3) oraz TRS (wariant 4). Widocznym jest, że zaledwie 


kilka próbek pochodzących z losowania Monte Carlo znajduje się w spodziewanym 


obszarze stanu granicznego. Można również zauważyć, że aproksymacja powierzchni w 


wariantach 2 1 3 przeprowadzona jest na podstawie punktów znajdujących się poza tym 


obszarem. 
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Rysunek 5.26 Schemat doboru punktów w czterech wariantach obliczeniowych 
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Rysunek 5.27 Schemat doboru punktów w czterech wariantach obliczeniowych 
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Jak wykazano w trakcie aproksymacji powierzchni odpowiedzi, dotaczenie 
próbek celowych do tych, które pochodzą z symulacji MC okazało się całkowicie zbędne. 
W dwóch pierwszych wariantach obliczeń aproksymacja powierzchni odpowiedzi 
cechowała się wysokimi błędami dopasowania, nieprzystającymi do wyników, jakie 
prezentowały pozostałe dwa warianty. Przyrosty całkowitego błędu aproksymacji (ERR) 
oraz błędu aproksymacji w przeliczeniu na liczbę użytych próbek (EPS) przedstawiono 
na rysunku 5.28. Wyniki uzyskane w przypadku zastosowania 1 i 2 wariantu 
obliczeniowego będą w niniejszym rozdziale podawane jedynie w celu wykazania 
błędów oszacowania, jakie pojawiają się w tych wariantach w stosunku do oszacowania 


metodą MC i techniką TRS. 


14 
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liczba próbek użytych do aproksymacji 


Rysunek 5.28 Całkowity błąd aproksymacji (ERR) oraz błąd aproksymacji w przeliczeniu na liczbę 
użytych próbek (EPS) wyników analizy metodą powierzchni odpowiedzi konstrukcji w czterech 
wariantach obliczeniowych 


W obu wariantach: 3 i 4, prezentujących zadowalająco mały poziom błędu 
oszacowania, interesującym zagadnieniem jest szybkość osiągania zbieżności wartości 
współczynników kierunkowych powierzchni odpowiedzi. Zmiany wartości 
współczynników kierunkowych, wyznaczonych na podstawie próbek symulacji Monte 
Carlo przedstawiono na rysunku 5.29. W przypadku wykorzystania 200 próbek 
pochodzących z symulacji Monte Carlo wartości współczynników kierunkowych 
stabilizują się po uwzględnieniu 55 próbek. Do tego momentu wahania uzyskanych 
rezultatów są niewielkie. Ponadto już przy minimalnej liczbie próbek (10-15) otrzymane 
wartości pozwalają przewidzieć rząd wielkości wartości współczynnika, jaki uzyska się 


po stabilizacji wyniku. 
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Rysunek 5.29 Stabilizacja wartości współczynników kierunkowych powierzchni odpowiedzi 
wyznaczonych na podstawie próbek symulacji Monte Carlo 


W przypadku aproksymacji powierzchni odpowiedzi na podstawie próbek 
wygenerowanych za pomocą techniki TRS (wariant 4) wartości współczynników 
kierunkowych zaczynają się stabilizować już po wykonaniu obliczeń 30 próbek 
(rysunek 5.29). Jednocześnie, do tego momentu, wahania wartości oszacowanych 
współczynników są stosunkowo duże i nie dają podstaw do określenia rzędu wielkości 


współczynnika tak szybko, jak było to możliwe w przypadku wariantu 3. 
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Rysunek 5.30 Stabilizacja wartości współczynników kierunkowych powierzchni odpowiedzi 
wyznaczonych na podstawie próbek wygenerowanych wg metody TRS 
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Wykorzystanie metody powierzchni odpowiedzi umożliwia dokonanie czytelnej 
graficznej interpretacji przebiegu zmienności mechanicznej odpowiedzi konstrukcji w 
zależności od wartości wejściowych zmiennych losowych zadania. Obrazowanie takie po 
pierwsze weryfikuje zasadność przyjęcia określonego stopnia wielomianu 
aproksymacyjnego, a po drugie dostarcza kluczowej wiedzy o topografii powierzchni 
odpowiedzi, jej ekstremach, a także o szacowanej wrażliwości odpowiedzi na wahania 
poszczególnych zmiennych wejściowych (można ją w przybliżeniu odczytać z nachyleń 
powierzchni). W przypadku dwóch zmiennych wejściowych obrazowanie takie jest 
możliwe w przestrzeni trójwymiarowej. 

Poniżej przedstawione zostaną graficzne postacie powierzchni odpowiedzi dwóch 
wariantów obliczeniowych o największej jakości wyników (wariantów 3 i 4). 

Na rysunku 5.31 przedstawiono powierzchnię odpowiedzi aproksymowaną w na 
podstawie 200 próbek uzyskanych metodą Monte Carlo (wariant 3). Współczynniki 
kierunkowe powierzchni wyznaczono przy użyciu programu RSM-Win. Po 
podstawieniu ich do wzoru 5.17 uzyskano następujące równanie aproksymowanej 


powierzchni: 


F(%,x,) =0,81976 + 3,5223x, +1,2775x, —4,4832x; —1,1537x> —2,1079x,x, (5.18) 


y(x; X) 


500 


x, [MN/m] 200 


100 
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Rysunek 5.31 Powierzchnia odpowiedzi uzyskana na podstawie 200 próbek losowania 
bezpośredniego Monte Carlo (wariant 3) 
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Z uwagi na lokalizację w przestrzeni zmiennych większości próbek 
wykorzystanych w aproksymacji (rys. 5.27) można oczekiwać, że przedstawiona na 
rys.5.31 powierzchnia najwierniej odwzorowuje zachowanie się konstrukcji w okolicy 


wartości średnich zmiennych k, i k,. Należy podkreślić, że w rozpatrywanym przypadku 


żadna z próbek wykorzystanych do aproksymacji powierzchni nie pochodziła z obszaru 
przekroczenia przez konstrukcję stanu granicznego a zaledwie jedna z nich znajdowała 
się w sąsiedztwie tego obszaru. Można zatem podejrzewać, że wyniki miar 
niezawodności uzyskane przy użyciu powierzchni aproksymowanej na podstawie próbek 
metody Monte Carlo nie są dokładne i należy poszukiwać w przestrzeni zmiennych 
próbek zapewniających dokładniejsze odwzorowanie powierzchni w obszarze stanu 
granicznego. 

W celu wyznaczenia przybliżonej powierzchni odpowiedzi w wariancie 4 
(rysunek 5.32) wykorzystano 200 próbek losowania ukierunkowanego. 
Po podstawieniu uzyskanych współczynników kierunkowych do wzoru 5.17 otrzymano 


następujące równanie powierzchni: 


$(x,,x,) = 0,34325+10,6953x, +1,7367x, —5,2915x; —6,3743x; —2,0104x,x, (5.19) 


200 
100 100 
X> [MN/m] x, [MN/m] 


Rysunek 5.32 Powierzchnia odpowiedzi uzyskana na podstawie 200 próbek losowania 
ukierunkowanego TRS (wariant 4) 


Należy podkreślić, że w odróżnieniu od wariantu wykorzystującego próbki 
metody Monte Carlo, przedstawiona powierzchnia nie obrazuje w sposób dokładny 


zachowania się konstrukcji w pełnym zakresie zmienności jej parametrów, a jedynie w 
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pobliżu funkcji stanu granicznego. Nachylenie powierzchni uzyskane tą techniką w 
znacząco lepszy sposób odwzorowuje obszar stanu granicznego zlokalizowany w okolicy 
„ogonów ” łącznego rozkładu zmiennych losowych zadania i pozwala na dokładniejsze 
odwzorowanie miar niezawodności rozważanego przykładu. 

Przedstawiony na rysunku 5.32 grzbiet powierzchni ma zasadność tylko dla 
obszaru po stronie mniejszych wartości realizacji zmiennych. Drugie zbocze grzbietu jest 


pomijalne - nie uczestniczy w obliczeniach miar bezpieczeństwa zagadnienia. 


Na rys.5.33 dodatkowo przedstawiono obie powierzchnie (uwidocznione osobno 
na rysunkach 5.31 oraz 5.32) ale tylko w zakresie wpływającym na estymację 
niezawodności konstrukcji. Charakter tych powierzchni jest wyraźnie zróżnicowany. Ze 
względu na liczbę punktów, na podstawie których obie powierzchnie są w tym obszarze 
aproksymowane, rozwiązanie uzyskane na podstawie metody TRS jest zdecydowanie 


bardziej wiarygodne. 


x, [MN/m] 150 


1 
100 00 120 140 160 


60 m x, [MN/m] 


Rysunek 5.33 Zestawienie fragmentów powierzchni odpowiedzi uzyskanych w wariantach: 3 i 4 w 
obszarze przecięcia z płaszczyzną y(x,,x,)=1. 


W dalszym etapie analizy, aproksymowane powierzchnie będą stanowiły podstawę do 


oszacowania miar bezpieczeństwa wieży. 


= 
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Analiza miar bezpieczeństwa wieży wykorzystująca metodę powierzchni 


odpowiedzi (RSM) 


Na podstawie aproksymowanych powierzchni odpowiedzi, w każdym z 4 
opisanych wcześniej wariantów doboru punktów próbkowania przeprowadzana została 
analiza niezawodności układu. W każdym z wariantów, dla kolejno zwiększającej się 
liczby próbek uwzględnionych w aproksymacji, wyznaczono prawdopodobieństwo 


awarii konstrukcji P, oraz wskaźnik niezawodności Hasofera — Linda (5,, ), 


f 
wykorzystujący gradienty współczynników kierunkowych powierzchni odpowiedzi 


(wzór (2.12)). Rezultaty przedstawionej analizy zaprezentowano na rys. 5.34 oraz 5.35. 
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Rysunek 5.34 Wyniki obliczeń wskaźnika niezawodności Hasofera — Linda ( /8,, ) na podstawie 
równań powierzchni odpowiedzi dla czterech wariantów zadania 
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Rysunek 5.35 Wyniki obliczeń prawdopodobieństwa awarii konstrukcji ( P, ) na podstawie równań 


powierzchni odpowiedzi dla czterech wariantów zadania 
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Zaprezentowane na rysunkach 5.34 oraz 5.35 rezultaty, poza oczywistym 
wnioskiem o konieczności odrzucenia wyników pochodzących z wariantów 1 i 2, 
wyraźnie wskazują, iż aproksymacja powierzchni odpowiedzi na podstawie próbek 
losowych metody Monte Carlo daje dobry (obarczony niewielkim błędem względnym) 
wynik miar bezpieczeństwa. Zbieżność tego wyniku jest jednak powolna - nie stabilizuje 
się on nawet po uwzględnieniu w aproksymacji liczby 200 próbek. Można również 
zaobserwować, iż rezultaty miar bezpieczeństwa są w tym przypadku wyższe niż w 
pozostałych wariantach. Wynikać to może z rozbieżności pomiędzy aproksymowaną i 
rzeczywistą odpowiedzią konstrukcji w przypadku małych wartości współczynników 
sprężystości podpór wieży kc oraz ky. 

Na tym tle, oszacowanie miar bezpieczeństwa konstrukcji wykorzystujące 
powierzchnię odpowiedzi aproksymowaną na podstawie próbek metody TRS gwarantuje 
wynik o znacznie lepszej jakości. Po pierwsze wartości /8,„, oraz P, obarczone są 
niewielkimi błędami względnymi oszacowania, mniejszymi niż w przypadku próbek MC. 
Ponadto rezultaty w tym przypadku szybciej uzyskują zbieżność - wyniki miar 
bezpieczeństwa można uznać za ustabilizowane po analizie 50 próbek. Należy również 
zwrócić uwagę , że powierzchnia odpowiedzi bardzo dobrze odwzorowuje rzeczywiste 
zachowanie się konstrukcji w obszarze stanu granicznego. Wynik jest zatem wiarygodny 
i nie charakteryzuje się zawyżeniem, cechującym rezultaty uzyskane w wariancie 3. 

Na podstawie przedstawionych rezultatów można wskazać, iż zaproponowana 
autorska procedura generacji próbek obliczeniowych w połączeniu z przyjętą metodą 
aproksymacji powierzchni odpowiedzi jest poprawna i dla przykładu rzeczywistej 
konstrukcji inżynierskiej osiąga lepsze rezultaty, niż techniki standardowe (Monte Carlo) 


lub popularne procedury związane z kryteriami odchyleń standardowych. 


5.3.3 Analiza wrażliwości konstrukcji na zmiany współczynników sprężystości 


podpór wieży kc, kr. 
Analiza wrażliwości konstrukcji na podstawie chmur rozproszenia wyników 


Pierwszym, najmniej czasochłonnym sposobem określenia wrażliwości układu na 
zmianę jego parametrów projektowych jest ilustracja obrazów rozproszenia rezultatów 
oraz określenia na ich podstawie wartości estymatorów współczynników korelacji 


liniowej (zgodnie z formułą 5.11). W przedstawionej analizie chmury obrazują wyniki 
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200 symulacji Monte Carlo. Zależności pomiędzy odpowiedzią konstrukcji a zmiennymi 
kc i ki przedstawiono na rysunkach: 5.36 oraz 5.37. Estymatory współczynników korelacji 
liniowej wyniosły odpowiednio: 7, „ = 0,8601 oraz 7, , = 0,423. 
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Rysunek 5.36 Kształt chmury rozproszenia wyników w przypadku zależności y(k.) 
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Rysunek 5.37 Kształt chmury rozproszenia wyników w przypadku zależności y(k, ) 


Rezultaty pobieżnej analizy wskazują, że zmiana parametru współczynnika 
sprężystości podpory pod słupem ściskanym (kc) w większym stopniu wpływa na 
odpowiedź konstrukcji —chmura wyników na rysunku 5.36 jest zauważalnie mniej 
rozproszona, korelacja pomiędzy punktami odpowiedzi jest zatem silniejsza niż ma to 


miejsce w przypadku chmury rezultatów przedstawionej na rys. 5.37. 
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Analiza wrażliwości konstrukcji -zastosowanie podejścia OAT (One-At-a-Time) 


Kolejnym sposobem określenia wrażliwości konstrukcji na zmiany jej 
parametrów projektowych — współczynników sprężystości podpór ke i k, było 
zastosowanie techniki OAT polegającej na przeprowadzaniu serii badań numerycznych 
uwzględniających zmienność tylko jednego z rozpatrywanych czynników. Aby tego 
dokonać, wygenerowano kolejne modele obliczeniowe wieży. Pierwsze 200 modeli 


wygenerowano zakładając, iż zmienna x=k 


c 


będzie zawsze równa 
k, = H, = 240 MN/m, natomiast zmienna x, =k, będzie losowana według rozkładu 
normalnego, którego wartość średnia i odchylenie standardowe wynoszą odpowiednio: 
Hy, =360 MN/m, o;, =100MN/m. Kolejne 200 modeli zakłada sytuację odwrotną — to 
zmienna x, =k, będzie zawsze stała i równa k, = 44, =360MN/m, natomiast zmienna 
x, =k, będzie przyjmowała wartości losowe o rozkładzie normalnym, którego wartość 


średnia oraz odchylenie standardowe wynoszą: My, = 240 MN/m , O, = 65 MN/m. 


W pierwszym z wymienionych przypadków obliczeniowych (stała kc) otrzymano 


następujące wartości ekstremalne mnoznika obciążenia dopuszczalnego: y,., =1,518 
oraz y,,,, =1,746. Rozstęp wyników (różnica pomiędzy maksymalnym i minimalnym 


rezultatem) wyniósł Ay,, =0,228. 


Na rysunkach 5.38 oraz 5.39 przedstawiono zmianę wartości średniej i odchylenia 
standardowego mnożnika dopuszczalnego obciążenia wieży w przypadku próbek analizy 


wrażliwości zmiennej x,=k,. Na rys.5.40 zamieszczono poglądowy histogram 


rezultatów mnożnika m,, uzyskanych dla rozważanych 200 próbek. 


Pobrano z mostwiedzy.pl 
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X1= const = 240 MN/m 
1,71 4 X2 - losowe 


w. średnia mnożnika dopuszczalnego 
obciążenia wieży Maop [-] 
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Rysunek 5.38 Proces zmiany wartości średniej mnożnika obciążenia dopuszczalnego wieży w 
przypadku gdy k, = 4, =240 MN/m = const 


xı = const = 240 MN/m 
X2 - losowe 


odch. stand. mnożnika dop.obciążenia 
wieży Mdop [-] 
© 
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Rysunek 5.39 Proces zmiany odchylenia standardowego mnoznika obciążenia dopuszczalnego 
wieży w przypadku gdy k, = i, = 240 MN/m = const 


60 PRZYPADEK: x,=k,=240 MN/m=const 


Hm... =1,685 [-] 


cop 


Gm, = 0.035 [-] 


dt 
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Rysunek 5.40 Histogram wartości mnożnika obciążenia dopuszczalnego wieży w przypadku gdy 
k, = Hg, = 240 MN/m = const 
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Obliczenia powtórzono dla przypadku, który zakłada stałą wartość współczynnika 
kı. Otrzymano następujące wartości ekstremalne mnożnika obciążenia dopuszczalnego: 


Jmin 71097 oraz Ving, 1821. Rozstęp wyników wyniósł Ay, =0,725, co jest 


wartością zbliżoną do rezultatu uzyskanego przy wykorzystaniu bezpośredniej metody 
Monte Carlo, w przypadku gdy obydwa współczynniki sprężystości traktowane są jako 
losowe. 

Narysunkach: 5.41 oraz 5.42 przedstawiono graficzną analizę zbieżności wartości 
średniej i odchylenia standardowego mnożnika obciążenia dopuszczalnego wieży, 
natomiast rys. 5.43 przedstawia histogram wyników mnożnika dopuszczalnego 
uzyskanych w rozważanym przykładzie. 
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Rysunek 5.41 Proces zmiany wartości średniej mnożnika obciążenia dopuszczalnego wieży w 
przypadku gdy k, = 4,, =360 MN/m = const 
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Rysunek 5.42 Proces zmiany odchylenia standardowego mnożnika obciążenia dopuszczalnego 
wieży w przypadku gdy k, = Hy, =360 MN/m = const 
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PRZYPADEK: x2=k,=360 MN/m=const 


Union = 1.667 [-] 
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Rysunek 5.43 Histogram wartości mnożnika obciążenia dopuszczalnego wieży w przypadku gdy 
k, = 44, =360 MN/m = const 


Na rys. 5.44 przedstawiono wykresy zależności pomiędzy wartościami 
współczynnika sprężystości podpory (podpór) wieży (założonego w danym wariancie 
obliczeniowym jako losowy) a dopuszczalnym mnożnikiem obciążenia wieży . Kolorem 
czerwonym zaznaczono 200 wartości mnożnika m, otrzymanych w przypadku gdy 


zmienna k, =, =240 MN/m = const., zaś kolorem zielonym oznaczono wyniki 


otrzymane dla 200 modeli, w których k, = 4, = 360 MN/m = const. 
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Rysunek 5.44 Zależności pomiędzy losową wartością współczynnika sprężystości podpory (podpór) 
wieży a zmienną odpowiedzi konstrukcji w dwóch wariantach obliczeniowych; linią przerywaną 
zaznaczono aproksymację wyników wielomianem stopnia czwartego 
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Przedstawione na wykresie zależności dostarczają informacji o zachowaniu się 
konstrukcji na skutek zmiany rozpatrywanego parametru. Na podstawie porównania 
szybkości zmiany wyników obliczeń oraz zakresu w jakim przyjmuje wartości 
odpowiedź konstrukcji w poszczególnym wariancie obliczeniowym można wnioskować 
o większej wrażliwości układu na zmianę parametru kc, co potwierdza wcześniejsze 
obserwacje. 

Na rysunku 5.45 przedstawiono te same zależności jednak w tym przypadku 
zmienne losowe zostały znormalizowane. Współczynnik nachylenia stycznych do 
aproksymowanych krzywych w punktach wartości średnich można interpretować jako 
kolejny wariant opisu wrażliwości .Uzyskane wykresy stanowiły podstawę obliczeń 


zawartych w dalszej części pracy. 
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Rysunek 5.45 Zależności pomiędzy znormalizowanymi zmiennymi losowymi z, i z, a zmienną 


odpowiedzi konstrukcji w dwóch wariantach obliczeniowych; linią przerywaną zaznaczono 
aproksymację wyników wielomianem stopnia czwartego 


Analiza wrażliwości wieży na podstawie różnic pomiędzy rozstępami wyników 


obliczeń 


Na podstawie rozstępów wyników uzyskanych w 3 wariantach obliczeniowych, 
opisanych w tabeli 5.11, sformułowano procedurę umożliwiającą oszacowanie 


współczynników wrażliwości obu zmiennych. 
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przypadek 


analizowany parametr k, - losowe k, =, = const k, -losowe 
k, -losowe k, = losowe k, = 44, = const 


minimalna wartość zmiennej 


0,83579 (-30) 


odpowiedzi - Vmin 


maksymalna wartość zmiennej 


odpowiedzi - Vinax 1,82175 (+30) 


Tablica 5.11 Ekstremalne wartości mnożnika obciążenia dopuszczalnego wieży oraz rozstępy 
wyników w trzech wariantach obliczeniowych uwzględniających losowość obydwu lub jednej ze 
zmiennych 


W przypadku współczynników wrażliwości dwóch zmiennych niezależnych zachodzi 
równość: 
m, +0," =1 (5.20) 
Podstawiając w miejsce obliczanych współczynników wrażliwości zmianę rozstepu 
pomiędzy obliczoną wartością przy zmiennej ustalonej, a wartością z pełnej analizy 
próbkowania bezpośredniego Monte Carlo, otrzymujemy: 
Ay,” + Ay,” =0,769” +0,265* = 0,6616 (5.21) 


Normalizując obie strony równania w ten sposób, aby spełnić równanie 5.20 uzyskamy: 


(n-Ay, ) +(n-Ay, | =(n-0,769) +(n-0,265) =n?-0,6616=1 (5.22) 


Zatem, n’ = =1,51152 oraz n=4/1,51152 =1, 22944, 


0,6616 
Podstawiając obliczoną wartość do wzoru 5.22 uzyskujemy następujące wartości 
wskaźników stanowiących miarę wrażliwości odpowiedzi konstrukcji na wahania 
zmiennych : 


a, =n- Ay, =1,229-0,769 =0,9454 ; a, =n-Ay, =1,229-0,265 =0,3258 (5.23) 


Analiza wrażliwości wieży na podstawie stycznych do wykresów wrażliwości w 


punkcie wartości średniej zmiennej odpowiedzi 


Kolejna propozycja oszacowania współczynników wrażliwości dla obydwu 
zmiennych zadania bazuje na wyznaczeniu stycznych do wykresów wrażliwości w 
punkcie wartości średniej zmiennej odpowiedzi. 


Styczne do wykresów krzywych, przedstawionych na rysunku 5.45, dane są równaniami: 
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ye” =0,1542-k, +1,96917 dla k, = m, =240 MN/m, k, -losowa (5.24) 

y” =0,0676-k, +1,96917 dla k, =, =360 MN/m, k, -losowa (5.25) 

Dalsza procedura postępowania jest identyczna jak w przypadku obliczeń 

przeprowadzonych na podstawie rozstępów wyników. Podstawiając w miejsce 

obliczanych współczynników wartości kątów nachylenia stycznych do wykresów 
wrażliwości w punkcie ich wartości średniej otrzymujemy: 

a + a, = 0,1542” +0,0676° = 0,0284 (5.26) 

Normalizując obie strony równania w ten sposób, aby spełnić równanie 5.20, uzyskamy: 


(n-a, ) +(n-a, | =(n-0,1542) +(n-0,0676) = n*-0,0284=1 (5.27) 


Zatem n = ./35,2766 = 5,93941. 


Podstawiając obliczoną wartość n do wzoru 5.27 otrzymujemy następujące wartości 
wskaźników, które można potraktować jako miarę wrażliwości mechanicznej odpowiedzi 
konstrukcji: 

a, =n: Ay, =5,939-0,154=0,916; a, =n-Ay, =5,939-0,0676= 0,4015 (5.28) 


Analiza wrażliwości wieży na podstawie krawędzi przecięcia powierzchni 


odpowiedzi z płaszczyznami x,y oraz x,y. 


Innym sposobem badania wrażliwości konstrukcji na wahania jej parametrów 
losowych była obserwacja zmiany wskaźnika niezawodności Hasofera — Linda ( B,,) 
wyznaczonego na podstawie równań powierzchni odpowiedzi. Aproksymacje 
przeprowadzono dla 100 wariantów modelu obliczeniowego, wygenerowanych według 
następującego schematu: 

= wybór zmiennej, która poddana zostanie analizie (np. x,=k,,) oraz ustalenie 
zbioru dyskretnych wartości, jakie przyjmować będzie analizowany czynnik. 
Założono, iż zmienna k,, równa będzie następującym wielkościom: 44 —30, ; 
My -2,75 O; > My 259 O; > My =2,29 O, My —20, > M 70%; My; My FOX; 
My, + 20, ; M, + 30, : 
Na podstawie wstępnej analizy zachowania konstrukcji ustalono, iż zakres 


( lt, —30; ; by, -20, ) obejmuje kluczowy obszar zmiany wskaźnika 
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niezawodności, dlatego też zdecydowano się na zagęszczenie wartości 
analizowanej zmiennej w tym przedziale; 

= dla każdej z realizacji zmiennej x, generacja z rozkładu normalnego pięciu wartości 
zmiennej x, =k,; 

« stworzenie modeli obliczeniowych wieży dla uzyskanych 50 par wartości 
zmiennych (x,,x, ); 

= powtórzenie procedury dla sytuacji odwrotnej — zmienna x, jest losowana z 
rozkładu normalnego, natomiast zmienna x, =k, przyjmuje dziesięć ustalonych 


odgórnie wartości. 


W przypadku każdej realizacji obliczono mnożnik obciążenia dopuszczalnego. Na 
podstawie wyników każdej z pięciopunktowych serii, wyznaczonych dla poszczególnych 
ustalonych wielkości analizowanej zmiennej, wykonano aproksymację krzywej 
odpowiedzi. W ten sposób opis ograniczono do jednej zmiennej losowej, a więc 
otrzymuje się pewnego rodzaju płaszczyznę przekrojową przez powierzchnię odpowiedzi 
zadania, co daje z kolei ślad przecięcia w postaci krzywej. 

Na podstawie równania w/w krzywej odpowiedzi, metodami gradientowymi oszacowano 


poszukiwany wskaźnik niezawodności /8,,. Dzięki odpowiednio ukierunkowanemu 
zestawowi płaszczyzn przekrojowych przez powierzchnię odpowiedzi możliwe jest 
graficzne przedstawienie zmiany wskaźnika niezawodności pomiędzy kolejnymi 
płaszczyznami a zatem kolejnymi dyskretnymi wartościami analizowanej zmiennej 
losowej. 

W pierwszej kolejności, zgodnie z wcześniej podanymi wytycznymi ustalono 
zbiór 10 dyskretnych wartości współczynnika sprężystości elementu podpierającego słup 
ściskany x, =k, =[45;61,25 ;77,5 ;93,75 ;110;175;240 ;305 ;370;435]MN/m, a 
następnie wylosowano 50 wartości zmiennej x, =k, o rozkładzie normalnym. W ten 
sposób uzyskano 50 realizacji wektora zmiennych x=[x,;x,] poprzez zestawienie 
każdej z dyskretnych wartości zmiennej x, z pięcioma losowo wygenerowanymi 
wartościami zmiennej x, . 

Na rys. 5.46 przedstawiono krawędzie przecięcia powierzchni odpowiedzi konstrukcji z 


10 płaszczyznami x,y zawierającymi kolejne ustalone wartości zmiennej x,. 
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Obok aproksymowanych krzywych odpowiedzi podano takze ich rownania dane w 


postaci wielomianów stopnia drugiego. 


1,90 - 
y=-205x2 +0,0018x> + 1,47 
1,80 4 vy = -106x+2 + 0,0013%x> + 1,51 
= -10-6X2 + 0,0013x + 1,45 
1,70 4 
y = -10:5x+2 + 0,0013x2 + 1,406 
1,60 4 = -106x22 + 0,0011x2 + 1,33 
1,50 4 


y=-10%x22 +0,0012x2 + 1,1677 
1,40 4 ss = y=-10x;2 + 0,0012x> + 1,11 
1,30 4 y =-20%x;2 + 0,0013x; + 1,032 
1,20 4 > BEE y =-10%5x;2 + 0,001x> + 0,97 
1,10 - 

oo —— Fo y = -607x22 + 0,0006x2 + 0,91 
1,00 4 


Xx1=|-30=45 MN/m O x,=H-2,750=61,25 MN/m 


o 
0,90 4 A x,=H-2,50=77,5 MN/m 
X 


mnożnik dopuszczalnego obciążenia wieży Mg [-] 


X  X1=H-2,250=93,75 MN/m 
0,80 4 x,=u-20=110 MN/m O X,=H|-0=175 MN/m 
0,70 4 x1=U=240 MN/m X x;,=uto=305 MN/m 
© X17U+20=375 MN/m —  X15U+30=435 MN/m 
0,60 7 r r 7 1 7 1 
80 180 280 380 480 580 680 780 


zmienna losowa - współczynnik spreżystości podpory x2=k;[MN/m] 


Rysunek 5.46 Kolejne przekroje przez powierzchnię odpowiedzi zagadnienia, wyznaczone dla 
ustalonych 10 dyskretnych wartości zmiennej losowej współczynnika sprężystości elementu 
podpierającego słup ściskany x, =k, 

Analogiczną procedurę przeprowadzono w przypadku gdy zmienna x, =k, 
przyjmuje wartości: x, =k, =[60;85 ;110 ;135;160 ;260;360 ;460 ;560 ;660] MN/m, 
natomiast zmienna x, losowana jest z rozkładu normalnego. Krawędzie przecięcia 
powierzchni odpowiedzi konstrukcji z 10 równoległymi płaszczyznami xy 
zawierającymi kolejne ustalone wartości zmiennej x, przedstawia rys. 5.47. Obok 


aproksymowanych krawędzi, na rysunku podano także równania krzywych odpowiedzi 


dane w postaci wielomianów stopnia drugiego. 
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ŚR y = -506x42 +0,0037x, + 1,16 
T igo] 6 y= -408x42 + 0,0032x4 + 1,16 
2 y = -60%x,2 +0,0045x, + 1,04 
£ y = -508x412 + 0,0037x; + 1,13 
z 1707 y=-30%xq2 + 0,0027x, + 1,1885 
g = -208x42 + 0,0022x4 + 1,03 
Z 160 4 y = -40-5x,2 + 0,003x4 + 1,07 
= y = -505x42 + 0,0034x, + 0,995 
3 1,50 4 = -50-6x2 +-0,0035x1-+-0,93 
$ 
o, 1,40 + y = -30$x12 + 0,0024x; + 0,914 
e 
© 1,30 4 
X 
2 1,20 - 
8 © X2="-30=60MN/m O x,=1-2,750=85MN/m 
a 110- A x>=h-2,50=110MN/m X x»=h-2,250=135MN/m 
N X x2=u-20=160MN/m O X,=p-0=260MN/m 
g 1,00 4 + x2=H=360MN/m - X,="+0=460MN/m 

— = x2=H+20=560MN/m + x,=u+30=660MN/m 

45 145 245 345 445 545 645 


zmienna losowa - współczynnik spreżystości x1=kc [MN/m] 


Rysunek 5.47 Kolejne przekroje przez powierzchnię odpowiedzi układu wyznaczone dla ustalonych 
10 dyskretnych wartości zmiennej losowej współczynnika sprężystości elementów podpierających 
słupy rozciągane x, =k, 

Wyniki zaprezentowane na rysunkach 5.46 oraz 5.47 zestawiono na wspólnym 
wykresie (rys. 5.48), dzięki czemu można zaobserwować, iż przedstawione 20 krzywych, 
wyrysowanych na podstawie 5 punktów obliczeniowych każda, są odzwierciedleniem 


powierzchni odpowiedzi konstrukcji. 


yX, X,) [-] 


ey D e eres a 


xk [MINA] ene x,=k, [MN/m] 


Rysunek 5.48 Wzajemna przestrzenna zależność zmiany krzywych odpowiedzi 
od ustalonej wartości zmiennych losowych x, i x, 
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Dla kazdej z 20 aproksymowanych krzywych odpowiedzi zaznaczonych na rysunkach 


5.46 - 5.47 wyznaczono wartość wskaźników niezawodności Hasofera — Linda £,,,, 


zgodnie z formułą 2.12. Wyniki szacowania wskaźników zestawiono w tablicy 5.12. 


x, - ustalone, x, - losowe 


x, - ustalone, x, - losowe 


x, [MN/m] Bu [-] 


x, [MN/m] 


Bn [-] 


Tablica 5.12 Wskaźniki niezawodności Hasofera — Linda ( B,, ) obliczone dla 20 przypadków 


krzywych odpowiedzi 


Na podstawie danych zebranych w tablicy 5.12, na rysunku 5.49 przedstawiono 


zależności pomiędzy wartościami współczynników sprężystości k,, k,, a wskaźnikiem 


niezawodności /5,, . 


74 6,728) _. 
ESP 

56,5 ee EF 6,767 
a 5,9447" 

6 - + 
g a 6,031 
555 | z 
1 y 
5 5o a 
g / © 4,823 
845 4 7 ° 4,517 5 > 
ET r et 
3 / ° 3,866-% 755 i i 
Q | i 
$35 T 3205 77 3,386 
G 37 are 3,281 © x, - state, xz - losowe 
(0) CA 
= 25 4 ‘a © x, - stałe, x, - losowe 
x 
$ 15 | ° 1,650 

; i 

i 
1 T T T T T T 1 
0 100 200 300 400 500 600 700 


zmienne: x4, X2[MN/m] 


Rysunek 5.49 Aproksymacja krzywych wrażliwości na podstawie wskaźników niezawodności 
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Hasofera — Linda ( 5,, ) 
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Na rysunku 5.49 linią przerywaną zaznaczono krzywe będące aproksymacją 
wyżej opisanych wartości wskaźników niezawodności dla odpowiednich zmian 
współczynników sprężystości podpór wieży k,, k,. Krzywe te będą stanowiły podstawę 


obliczeń wrażliwości, przeprowadzonych w dalszej części pracy. 


Analiza wrażliwości niezawodności konstrukcji na podstawie stycznych do 


wykresów wrażliwości w punkcie wartości średniej zmiennej odpowiedzi 


Styczne do wykresów krzywych wrażliwości przedstawionych na rys. 5.49 
opisują następujące równania: 
ye” =0,460-k, +6,031 dla k, = w, = 240 MN/m, k,- losowe (5.29) 
y =0,117-k, +4,415 dla k, = u, =360 MN/m, k,- losowe (5.30) 
Podstawiając w miejsce obliczanych wskaźników wrażliwości równania 5.20 wartości 
współczynników kierunkowych stycznych do wykresów wrażliwości w punkcie wartości 
średniej zmiennej odpowiedzi, otrzymujemy: 
a,” +a,° =0,460° +0,117° =0,225289 (5.31) 
Normalizując obie strony równania tak aby spełnić równanie 5.20 uzyskamy: 
(n-a, ) +(n-a, } =(n-0,460) +(n-0,117) =n*-0,225289=1 (5.32) 
Zatem n =„./4,438743 =2,10683. 


Podstawiając obliczoną wartość n do równania 5.32 otrzymujemy następujące wartości 


a, oraz Ot, : 


a, 


kę 


=n: Ay, =2,1068-0,46=0,9691; a, =n-Ay, =2,1068-0,117=0,2465 (5.33) 


Wielkości te można traktować jako przybliżenie wskaźników wrażliwości niezawodności 


konstrukcji. 


Analiza wrażliwości odpowiedzi konstrukcji- wyniki uzyskane w programie 


komercyjnym COMREL 


Kolejno, problem zaimplementowano w środowisku komercyjnego programu do 
obliczeń niezawodnościowych i probabilistycznych StruRel— COMREL 9.0 RCP GmbH. 


Do pamięci programu wprowadzono funkcję stanu granicznego daną równaniem 5.14 a 
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także wszystkie niezbędne dane probabilistyczne zmiennych k, i k,. Do obliczeń 


współczynników wrażliwości posłużono się techniką Hasofera — Linda. Wyniki działania 
programu, przedstawione w formie diagramu kołowego (domyślny sposób wyświetlania 


współczynników wrażliwości w programie), zaprezentowano na rysunku 5.50. 


kt 10.30 
KC 0.95 
Sum of a? 1.00 


Rysunek 5.50 Analiza wrażliwości wykonana w środowisku komercyjnego programu 
StruRel — COMREL 9.0 RCP GmbH 


5.3.4 Porównanie rezultatów analizy wrażliwości 


W tabeli 5.13 zestawiono wartości wskaźników wrażliwości odpowiedzi 
konstrukcji czterech niezależnych procedur obliczeniowych. W przypadkach procedur 
polegających na badaniu różnic między rozstępami wyników czy też nachyleniami 
stycznych dodatkowo podano różnice względne pomiędzy uzyskanymi w nich 
rezultatami a wartościami otrzymanymi w programie komercyjnym. 


wskaźniki wrażliwości wyznaczone na podstawie: 


różnic w różnic w 
różnic w nachyleniu nachyleniu 
rozstępach stycznych stycznych 
wyników do wykresów do wykresów 
mnożników m wskaźników 8 


zmienna losowa programu 


COMREL 


dop 


0,9159 0,9691 

A=3,59% A=2,01% 

0,4015 0, 2465 

A=33,83% A=17,83% 

Tablica 5.13 Porównanie wartości współczynników wrażliwości zmiennej odpowiedzi - mnożnika 
obciążenia dopuszczalnego wieży 
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Podsumowując, można stwierdzić, że wrażliwość odpowiedzi konstrukcji 
(wyrażonej zarówno poprzez mnożnik obciążenia dopuszczalnego, jak również poprzez 
wskaźnik niezawodności) na zmianę współczynnika sprężystości podpory słupa 
ściskanego jest wyraźnie większa, niż na zmianę współczynników sprężystości podpór 
słupów rozciąganych. 

Ponadto, wyniki zaprezentowane w tablicy 5.13 uzasadniają stwierdzenie, iż 
autorskie propozycje oceny wrażliwości są poprawne. Uzyskane rezultaty w sposób 
właściwy wskazują parametr, którego wpływ na mechaniczne zachowanie się układu jest 
większy. Dodatkowo rząd uzyskanych wielkości pozostaje ten sam a błędy oszacowania 
wrażliwości na zmianę zmiennej kluczowej w problemie względem wyników 
uzyskanych w programie komercyjnym są niewielkie (nie przekraczają 5%). 

Badanie wrażliwości układu na drodze analizy graficznej reprezentacji 
zachowania się konstrukcji prowadzi do właściwego wniosku, na wahania którego z 
parametrów zmienna wyjściowa jest wrażliwsza. Należy jednak podkreślić, że liczba 
realizacji, na podstawie której dokonuje się obrazu odpowiedzi konstrukcji musi być 
wystarczająca (jest bowiem uzależniona od stopnia skomplikowania obserwowanej 
odpowiedzi, tj. stopnia wielomianu aproksymującego) oraz dobrana w sposób świadomy 
i gwarantujący objęcie analizą całego krytycznego obszaru. Zaproponowana procedura 
autorska wykorzystująca punkty obliczeniowe w pobliżu funkcji stanu granicznego 


wydaje się spełniać te wymagania. 
5.4 Kopuła — kratownica przestrzenna 


Kolejnym przykładem wykorzystania metod probabilistycznych jest analiza 
wpływu imperfekcji geometrycznych, wynikających z obowiązujących w praktyce 
inżynierskiej tolerancji w zakresie niedokładności wykonawczych. Jako przykład 
wybrano małowyniosłą kopułę kratownicową, której stany graniczne są niezwykle 
wrażliwe na zmianę poszczególnych parametrów modelu. Większość umieszczonych w 


rozdziale obliczeń i wniosków prezentowano m.in. w [Sorn i inni, 2015, 2016]. 
5.4.1 Przedstawienie modelu 


Analizowana konstrukcja kratownicy przestrzennej, opisana w pracy [Rakowski i 
inni, 1993], składa się z 24 elementów o przekroju rurowym , wykonanych ze stali S355. 


Całkowita średnica konstrukcji wynosi 50 m. Najwyższy jej węzeł (nr 13) położony jest 


Pobrano z mostwiedzy.pl 
J 
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na wysokości 8,216 m ponad poziomem podparcia. Elementy kratownicy połączone są 
przegubami kulowymi. Geometrię analizowanej konstrukcji przedstawiono na rysunku 
5.51. Obliczenia wykonano za pomocą programu MSC NASTRAN przy założeniu 


geometrycznej i materiałowej nieliniowości. 


Rysunek 5.51 Przestrzenna konstrukcja kratownicowa- geometria [Sorn i inni, 2015] 
5.4.2 Wybór metody obliczeniowej 


W pierwszej kolejności wykonano obliczenia wstępne konstrukcji, mające na celu 
testowanie kilku różnych metod probabilistycznych i wyłonienie spośród nich procedury 
optymalnej. W tym celu założono prosty przypadek obciążenia jakim jest pojedyncza, 
pionowa siła skupiona przyłożona w węźle nr 13. Przyjęte obciążenie, mimo iż nie jest 
realistyczne, ma duży wpływ na odpowiedź mechaniczną układu i umożliwia obserwację 
zachowania modelu na skutek zmiany założonych parametrów projektowych. 
Pierwszym krokiem było wyznaczenie wartości obciążenia niszczącego, które wyniosło 


R=2612,3kN. Obliczenia probabilistyczne wykonano przy założeniu, że konstrukcję 


opisują dwie zmienne losowe: geometryczna - przemieszczenie u wierzchołka 
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kratownicy (węzeł 13 na rys. 5.51) oraz materiałowa — moduł Young’a E. Przyjęto 
następujący opis imperfekcji geometrycznych: wartość średnia m, =0,00m oraz 
odchylenie standardowe o, =0,08m . Parametry rozkładu normalnego zostały tak 
dobrane, aby zakres wygenerowanych przemieszczeń węzła nr. 13 mieścił się w 
granicach (—0,13+0,13)m. Przemieszczenia w podanym zakresie odpowiadają 
imperfekcjom wywołanym przez zmiany długości prętów kratownicy. Wahania 
wymiarów elementów (tolerancje) są następstwem niedokładności pojawiających się w 
procesie produkcyjnym i są regulowane przez normy. Początkowe przemieszczenia 
węzłów 7 — 12 zostały przyjęte proporcjonalnie do wygenerowanych imperfekcji węzła 
13-tego. Moduł Young’a został także opisany rozkładem Gaussa, ze średnią 
m, =210,0GPa oraz odchyleniem standardowym o, =4,0GPa. Przykładową 


generację zmiennych u i E przedstawiono na rys. 5.52 
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Rysunek 5.52 Generacja dwóch zmiennych losowych — imerfekcji geometrycznych i modułu 
Young'a 


m= 


Obliczenia wykonano za pomocą kilku metod. W przypadku bezpośredniej 
metody Monte Carlo w wygenerowano 60 realizacji i na tej podstawie uzyskano opis 


obciążenia niszczącego: jego wartość średnią m, = 2653,0kN , odchylenie standardowe 
0,=232,54kN oraz skośność 7,=0,1004kN. Wyniki te traktowane są jako 


rozwiązanie dokładne, pomimo tego, że 60 realizacji może wydawać się próbą 
niewystarczającą. Dla porównania rozrzutów wyników dodatkowo wygenerowano 36, 


24, 16 oraz 9 par zmiennych losowych. Wyniki umieszczono w tablicy 5.14. 
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Liczba 
Zastosowana metoda realizacji 


Odchylenie standardowe 
o, [kN] 


wartość średnia 
m; [kN] 


Monte Carlo 2653,21 232,54 
Monte Carlo 2658,71 229,21 
Próbkowanie warstwowe — równe 

przedziały 6x6 

Próbkowanie warstwowe — równe 

prawdopodobieństwa 6x6 

Monte Carlo 2651,48 227,16 
Próbkowanie warstwowe — równe 
przedziały 5x5 

Próbkowanie warstwowe — równe 
prawdopodobieństwa 5x5 

Monte Carlo 2602,15 253,78 
Próbkowanie warstwowe — równe 
przedziały 4x4 

Próbkowanie warstwowe — równe 
prawdopodobieństwa 4x4 

Monte Carlo 2728,93 309,69 
Próbkowanie warstwowe — równe 
przedziały 3x3 

Próbkowanie warstwowe — równe 
prawdopodobieństwa 3x3 

PEM 2620,48 229,38 


2637,43 388,28 


2604,47 255,01 


2645,16 410,19 


2630,21 237,22 


2647,99 389,89 


2614,67 299,57 


2633,06 348,67 


2656,06 328,89 


Tablica 5.14 Obciążenie graniczne kopuły kratownicowej — analiza z uwzględnieniem dwóch 
zmiennych losowych: przemieszczenie węzła nr 13 oraz modułu Young’a [Sorn i inni, 2015] 


Zastosowanie metody Monte Carlo jest niezwykle pracochłonne dlatego też 
poszukiwane są możliwości ograniczenia liczby realizacji poprzez zastosowanie technik 
redukcyjnych. Pierwszą z zastosowanych metod było próbkowanie warstwowe z liczbą 
podziałów przestrzeni zmiennych równą kolejno: 6x6, 5x5, 4x4 1 3x3. W ten sposób 
nawiązano do obliczeń wykonanych metodą Monte Carlo z wykorzystaniem tej samej 
liczby próbek. Obliczenia metodą próbkowania warstwowego wykonano w dwóch 
wariantach — pierwszy z nich przewidywał podział przestrzeni na warstwy o takich 
samych wymiarach, drugi zaś na warstwy o równych prawdopodobieństwach znalezienia 
się w danej podprzestrzeni. Wyniki umieszczono w tablicy 5.14. Pobieżna analiza 
uzyskanych rezultatów wskazuje, że te same wielkości wyznaczone innymi metodami 
różnią się od siebie w sposób znaczący, przez co trudno zdecydować, która z 


zastosowanych metod jest optymalna oraz jaką liczbę realizacji należy uznać za 
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minimalną. Co więcej w przypadku gdy przyjęta jest większa liczba zmiennych losowych 
wszystkie te metody będą mało efektywne. Z tego powodu kolejne obliczenia wykonano 
za pomocą metody PEM (tablica 5.15), która w znacznym stopniu ogranicza liczbę 
realizacji — w przypadku dwóch zmiennych losowych obliczenia należy wykonać jedynie 
dla czterech modeli deterministycznych. Wyniki uzyskane za pomocą tej metody 
wykazują niewielki błąd w stosunku do poziomu odniesienia, a mianowicie 1,2% w 
stosunku do wartości oczekiwanej i 1.4% w przypadku odchylenia standardowego. Ze 
względu na niewielkie różnice wyników oraz znaczną redukcję liczby realizacji w 


dalszych obliczeniach wykorzystano wyłącznie metodę PEM. 


JPG) 
697,545 1946276 
610,8075 1492343 
724,575 2100036 


587,548 1380848 


E[R]=2620,475kN, Var[R] = 6919503 — 2620,475° =52614,07kN”,o, = 229,38 kN 


Tablica 5.15 Rezultaty obliczeń w przypadku dwóch zmiennych losowych — metoda estymacji 
punktowej PEM [Sorn i inni, 2015] 


5.4.3 Analiza wrażliwości 


Dalszą analizę wykonano dla realistycznego obciążenia kratownicy śniegiem. Na 
podstawie normy PN-EN 1990-1-3 przyjęto następujące wzory opisujące obciążenie : 
s=LCCs, [KN/m?], (5.34) 
gdzie: 
s, — wartość charakterystyczna obciążenia śniegiem gruntu w rozpatrywanym miejscu 
[kN/m], 


H, — współczynnik kształtu dachu, 
C, współczynnik ekspozycji (przyjęto C, =1,0), 
C 


, —współczynnik termiczny (przyjęto C, =1,0). 
Rozpatrzono dwa przypadki — symetrycznego i niesymetrycznego obciążenie śniegiem 


(rys.5.53). 
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Rysunek 5.53 Symetryczne i niesymetryczne obciążenie śniegiem [Sorn i inni, 2015] 


W opisie kratownicy zastosowano cztery zmienne losowe: początkowe pionowe 


przemieszczenie węzła 13 (x,), przemieszczenie węzłów 7-14 (x,), Moduł Young’a 
(x;) oraz pole powierzchni przekroju poprzecznego prętów (x,). Przemieszczenia 


węzłów 7-14 opisano w identyczny sposób jak węzła 13. Dane dotyczące modułu 
Young 'a także pozostawiono takie same jak w analizie dwóch zmiennych losowych. 
Średnią oczekiwaną przekroju poprzecznego oraz jego odchylenia standardowe opisano 


rozkładem normalnym o następujących parametrach: m, =394,458 cm? oraz 


On, =2,0cm*. Zastosowano wersje Honga metody estymacji punktowej (PEM), 


omówioną w rozdziale 3.3 pracy. Warto podkreślić, że wersja Honga metody PEM, 
stosowana w przypadku przyjętej dużej liczby zmiennych losowych, może dawać błędne 
wyniki [Christian i inni, 2002], jednak dla 4 zmiennych jej stosowanie wydaje się 
uzasadnione. Przebieg obliczeń wykonanych dla symetrycznego i niesymetrycznego 
obciążenia śniegiem przedstawiono w tabelach 5.16 oraz 5.17 [Sorn i inni, 2015]. Należy 
podkreślić, że w analizie każdego z wariantów obciążenia wystarczające było 


przeprowadzenie obliczeń dla ośmiu deterministycznych modeli kratownicy. 
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y 
[kN /m?] 


X3 mean X4 mean 
=210,0 | =394,458 


X2, mean X3, mean X4, mean 


=0,16 | =0,0 | =210,0 |=394,458 


,mean X X3, mean X4, mean 


0,0 |=-0,06 | =210,0 | =394,458 


,mean Z X3, mean X4, mean 


0,0 =0,06 =210,0 | =394,458 
x ,mean X2, mean — X4 mean 


= 0,0 = 0,0 = 394,458 


X ,mean X2 mean X4 mean 
=0,0 = 0,0 = 394,458 
X4 


X ,mean X2 „mean 


=0,0 =0,0 = 354,458 


x „mean Xz, mean X44 


=0,0 | =0,0 = 434, 458 


Tablica 5.16 Symetryczne obciążenie śniegiem - rezultaty obliczeń w przypadku 4 zmiennych 
losowych — metoda estymacji punktowej Honga [Sorn i inni, 2015] 


X mean X3 „mean X4, mean 


=0,0 | =210,0 | =394,458 


X2, mean X3 „mean X4, mean 
=0,0 | =210,0 | =394,458 


X, mean A. x3. mean X4, mean 


=0,0 | =-0,06 | =210,0 | =394,458 


x, mean Xay X3, mean X4, mean 


=0,0 | =0,06 | =210,0 | =394,458 


AT, mean Xa, mean X3 X4, mean 


X, mean Xz, mean X34 X4, mean 


=0,0 | =0,0 | =218,0 | =394,458 


x), mean Xa, mean X3, mean ZIE 


=0,0 | =0,0 |=210,0 | =354,458 


X, mean X, mean X3 „mean X44 


= 0,0 =0,0 | =210,0 | =434,458 


4,6039 


E[R]=4,6039kN, Var|R] = 21,642 —4,6039* = 0,446kN* „o, = 0,668kN 


Tablica 5.17 Asymetryczne obciążenie śniegiem — przypadek 4 zmiennych losowych — metoda 
estymacji punktowej Honga [Sorn i inni, 2015] 
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Na podstawie obliczen uzyskano nastepujace wartosci oczekiwane i odchylenia 
standardowe obciążenia granicznego Sniegiem, dla obciążenia symetrycznego — 
23,106kN/ m° i 5,43kN/m” oraz niesymetrycznego — 4,6039kN/m” i 0,668kN/m”. 
Na szczególną uwagę zasługuje niezwykle mała wartość oczekiwana i odchylenia 
standardowego w przypadku niesymetrycznego obciążenia śniegiem. 

Estymacja niezawodności kratownicy wymaga przyjęcia odpowiedniego rozkładu 
obciążenia. Na podstawie wieloletnich pomiarów, przedstawionych np. w [Gwóźdź i 
Machowski, 2011] zastosowano rozkład Gumbela jako optymalny dla tego typu opisu 
[Sorn i inni, 2016]: 

f (x) =4,95:exp|-4,95(x—0,259) ]-exp| —exp(—4,95(x-0,259))| (5.35) 


Postać histogramu oraz aproksymowanej na jego podstawie funkcji gęstości rozkładu 
prawdopodobieństwa przedstawiono na rys. 5.54 


3 


N 


częstość [-] 


— 
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Rysunek 5.54 Histogram oraz przybliżona funkcja gęstości prawdopodobieństwa obciążenia 
śniegiem (rozkład Gumbela o parametrach: 0, = 4,95, u=0,259) 


W celu opisu wrażliwości niezawodności konstrukcji na zmianę przekroju 
poprzecznego elementów kopuły wykonano szereg obliczeń, których wyniki w 
przypadku symetrycznego obciążenia śniegiem, przedstawiono na rys. 5.55. Na tej 


podstawie wyznaczono wskaźnik niezawodności modelu stosując wzór Hasofera-Linda. 


N/N MOST 


126 Rozdziat 5 Probabilistyczna analiza wrazliwosci - przyktady numeryczne 


2 T T T T m = "=P" T T 


rozkład sym. obciążenia śniegiem 


——— rozkład obciążenia granicznego 


A=100cm? 


A=150¢m? 


częstość [-] 


-1 3 7 11 15 19 
obciążenie / obciążenie graniczne 


Rysunek 5.55 Rozkłady: symetrycznego obciążenia śniegiem oraz odpowiedzi układu w zależności 
od przyjętego pola poprzecznego A, elementów 


Aby ułatwić opis zmienności niezawodności modelu wyniki zostały w 
standardowy sposób aproksymowane następującymi funkcjami [Sorn i inni, 2016]: 
— w przypadku obciążenia symetrycznego śniegiem (otrzymano współczynnik korelacji 
0,9975): 
p, =4.2965-10" A7 (5.36) 
— w przypadku obciążenia niesymetrycznego Sniegiem (współczynnik korelacji wyniósł 
0,99): 
Py = 45319-10743" (5.37) 


Wyniki graficznie przedstawiono na rys. 5.56. 
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Rysunek 5.56 Aproksymacja wyników obliczeń niezawodności kopuły[Sorn i inni, 2016] 


Zgodnie z Normą PN-EN-1990 poziom niezawodności jest Ściśle związany z 
przyjętym okresem ,,zycia” konstrukcji. Zakładając ten okres na poziomie 50 lat i klasę 


RC3 niezawodności, index 8,, powinien wynosić 4,3. Analizując rys. 5.56 łatwo można 


stwierdzić, że w przypadku obciążenia niesymetrycznego odpowiada to przekrojom 
poprzecznym A=300 cm”, podczas gdy dla obciążenia symetrycznego wystarczy poziom 
około A=100 cm’. 

Przedstawioną analizę można potraktować jako algorytm postepowania w 
określeniu wrażliwości niezawodności konstrukcji na zmianę pojedynczego parametru, 


zakończony konkretnymi wnioskami projektowymi. 
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W pracy przedstawiono probabilistyczną analizę wrażliwości mechanicznej 
odpowiedzi konstrukcji z uwagi na odchyłki geometryczne i zmianę własności 
materiałów. Skoncentrowano się na zaprezentowaniu skutecznych metod, łatwych w 
zastosowaniach inżynierskich i niewymagających zaawansowanej wiedzy lub 
specjalistycznego oprogramowania. 

Głównym elementem analizy jest estymacja niezawodności konstrukcji i na tej 
podstawie szacowanie jej wrażliwości. Budując algorytmy obliczeniowe przyjęto 
podstawowe założenie, że funkcja stanu granicznego nie jest znana i należy ją wyznaczyć 
na drodze iteracyjnej. Taka sytuacja ma miejsce w przypadku większości nawet 
najbardziej standardowych konstrukcji. Zagadnienia, w których funkcja stanu 
granicznego jest możliwa do zdefiniowania na drodze analitycznej, dotyczą wyłącznie 
pojedynczych elementów konstrukcji i ograniczonej klasy zagadnień. 

W pracy rozwinięto powszechnie znane i stosowane algorytmy obliczeń: metodę 
Monte Carlo (MC) wraz z metodami redukcyjnymi próbkowania warstwowego (SS) i 
hipersześcianu łacińskiego (LHS), metodę estymacji punktowej (PEM) oraz metodę 
powierzchni odpowiedzi (RSM). Obliczenia przeprowadzone w pracy wyróżnia 
wszechstronność analizy porównawczej i próba zdefiniowania metody optymalnej. 

W przypadku metody powierzchni odpowiedzi (RSM) skoncentrowano się na 
sposobie dobierania próbek, tak aby skutecznie wymodelować powierzchnię stanu 
granicznego. W tym celu przeprowadzono obliczenia dla wielu wariantów symulacji, 
także niestandardowych, np. łącząc próbki wykorzystane w metodzie estymacji 
punktowej (PEM), z próbkowaniem metody Monte Carlo (MC). Dokonano 
szczegółowego porównania wyników tych obliczeń. Zadawalające rezultaty uzyskano 
przy zastosowaniu metody próbkowania warstwowego (SS), odpowiednio klasyfikując 
próbki wygenerowane za pomocą metody MC. Wykazano, że dobór metody 
obliczeniowej zależy od analizowanego problemu, a także od spodziewanego poziomu 
niezawodności. Przy inżyniersko prawidłowo zdefiniowanym systemie zawodność jest 
określana liczbą na poziomie 10% i jej estymacja jest niezwykle trudna. Gdy konstrukcje 
charakteryzują się większa zawodnością obliczenia są zdecydowanie łatwiejsze. 
Określenie tego poziomu i na tej podstawie obranie właściwej metody estymacji 


niezawodności powinno być pierwszym krokiem wykonywanej analizy. 
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W celu optymalnego wykorzystania metody powierzchni odpowiedzi (RSM) w 
środowisku MATLAB stworzono program pozwalający na poszukiwanie próbek w 
otoczeniu stanu granicznego konstrukcji. Autorski algorytm próbkowania 
ukierunkowanego (TRS) jest najważniejszym elementem pracy. Technika ta umożliwia 
redukcję nakładu obliczeń, gdyż powierzchnia stanu granicznego jest estymowana tylko 
w pobliży funkcji granicznej. Takie podejście prowadzi do szybkiego ustabilizowania 
wyników. W pracy jedynie zainicjowano analizę algorytmu TRS. Konieczne są dalsze 
prace nad odpowiednią doborem próbek, np. identyfikacją i odrzuceniem tych, które 
zaburzają budowę wielowymiarowej funkcji. 

Najważniejszym elementem pracy jest analiza wrażliwości mechanicznej 
odpowiedzi konstrukcji na zmianę parametrów materiałowych oraz odchyłek 
geometrycznych. W pracy przetestowano kilka wariantów obliczeń. W pierwszym 
przykładzie — kratownicy von Misesa — z uwagi na teoretyczny model kratownicy 
uzyskano histogramy nośności granicznej w funkcji poszczególnych parametrów 
definiujących tę konstrukcję. Na tej podstawie określono kilka globalnych (dotyczących 
całego histogramu) oraz lokalnych (w punktach projektowych) współczynników 
wrażliwości. Parametry globalne (wartość oczekiwana i odchylenie standardowe) w 
sposób ogólny charakteryzują histogram stanu granicznego. Parametry lokalne, 
definiowane na podstawie krzywej aproksymowanej w otoczeniu punktu projektowego 
pozwalają na liczbową charakterystykę wrażliwości konstrukcji i umożliwiają 
porównanie wpływu poszczególnych zmiennych opisujących model. W celach 
porównawczych wyznaczono także wskaźniki wrażliwości Sobola. 

Należy podkreślić, że podejście, którego podstawę stanowi wynikowy histogram 
stanu granicznego nie jest możliwe w przypadkach bardziej skomplikowanych 
konstrukcji, wymagających znacznie większego czasu obliczeniowego dla pojedynczej 
wygenerowanej próbki. Zaproponowano więc rozwiązania pośrednie wykonując analizy 
parametryczne i badając wpływ zmiany poszczególnych parametrów na mechaniczna 
odpowiedź modelu. 

Kompleksowe obliczenia dla drugiego teoretycznego przykładu — kolumny 
Ziglera — pozwoliły na wstępne określenie wrażliwości parametrów poprzez analizę 
chmury rozproszenia wyników obliczeń w funkcji stałych opisujących elementy 
sprężystych połączeń modelu. W celu udokładnienia tych obliczeń zastosowano metodę 


powierzchni odpowiedzi przy wykorzystaniu próbek Monte Carlo w wersji bezpośredniej 
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oraz stosując kilka metod redukcyjnych. Najbardziej skuteczna okazała się metoda 
próbkowania warstwowego (SS). Wykazano, że równie skuteczne obliczeniowo jest 
autorska metod próbkowania ukierunkowanego (TRS). 

Najbardziej rozbudowanym przykładem jest analiza rzeczywistej wieży 
telekomunikacyjnej. W obliczeniach wykorzystano wcześniejsze doświadczenia, 
dotyczące zarówno szacowania wpływu poszczególnych parametrów na mechaniczną 
odpowiedź konstrukcji jak i charakteru tego wpływu. 200 realizacji wykorzystanych w 
metodzie Monte Carlo (MC) nie pozwoliło na określenie niezawodności konstrukcji. 
Cztery warianty obliczeń wykonane za pomocą metody powierzchni odpowiedzi (RSM) 
udowodniły, że najbardziej optymalną jest metoda próbkowania ukierunkowanego 
(TRS). Wyniki tej kompleksowej analizy pozwoliły na określenie wrażliwości 
konstrukcji z uwagi na zmianę parametrów jej podparcia. Obszerne wyniki wykorzystano 
w aproksymacji funkcji opisujących wrażliwość, a następnie wyznaczenie odpowiednich 
liczbowych parametrów. Przyjęto autorski opis współczynników wrażliwości bliski 
wynikom analizy wykonanej za pomocą komercyjnego programu. Końcowym efektem 
obliczeń jest powierzchnia wrażliwości, którą jednak można uzyskać jedynie w wyniku 
czasochłonnych obliczeń w obszarze bliskim punktu projektowego. Jej znaczenie jest 
więc ograniczone. 

W ostatnim przykładzie zaprezentowano kompleksowe obliczenia kopuły 
prętowej. Stwierdzono, że w tym przypadku optymalną metodą jest metoda estymacji 
punktowej (PEM). Analiza wrażliwości pozwoliła na określenie minimalnego pola 
przekroju poprzecznego prętów przy założeniu 50-cio letniego normowego poziomu 


niezawodności. 
Reasumując, najważniejszymi elementami pracy są: 


= prezentacja metod wstępnej charakterystyki problemu, umożliwiających 
porównanie wpływu poszczególnych parametrów na mechaniczną odpowiedź 
konstrukcji oraz określenie charakteru tego wpływu, 

u wszechstronne przetestowanie metody powierzchni odpowiedzi (RSM) z uwagi 
na dobór próbek, pozwalających na właściwą aproksymację funkcji stanu 


granicznego, 
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= zaproponowanie autorskiej metody ukierunkowanego doboru próbek (TRS), 
który wydaje się być najwłaściwszą metodą estymacji niezawodności dowolnego 
wielowymiarowego problemu losowego, 


= zdefiniowanie opisu wskaźników wrażliwości konstrukcji. 
Kierunkiem dalszego rozwoju proponowanych metod jest: 


= udokładnienie metody próbkowania ukierunkowanego z eliminacją punktów 
zaburzających algorytm estymacji powierzchni stanu granicznego, 

« klasyfikacja metod obliczeniowych w zależności od cech charakterystycznych 
problemu, 


= ujednolicenie parametrów (wskaźników) opisujących wrażliwość problemu. 
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